IN8008 Organisatorisches und erste Wintersemester 2013 ... · PDF fileIN8008, Wintersemester...

Scientific Computing in Computer Science, Technische Universitat Munchen

Einfuhrung in die wissenschaftlicheProgrammierung

IN8008

Tobias Neckel

Wintersemester 2013/2014

T. Neckel: Einfuhrung in die wissenschaftliche Programmierung

IN8008, Wintersemester 2013/2014 1


Teil I

Organisatorisches und ersteSchritte




Inhalte und ZieleEinfuhrung in die• Keine Programmiererfahrung notig• Umgang mit PC/Betriebssystem sollte bekannt sein• Mischung aus einfachen und schwierigen Teilen• Kratzt an vielen Stellen nur an der Oberflache

wissenschaftliche• Mathematisches und physikalisches Grundverstandnis• Beispiele aus dem wissenschaftlichen Bereich• Hauptziel bleibt Verstandnis der Programmierkonzepte

Programmierung• Hauptsachlich Implementierung• Aber auch der Entwurf von Programmen• Ubliche Programmierkonstrukte und wissenschaftliche

Bibliotheken• Qualitat




Organisatorisches

Raum und Zeit• Modul IN8008• Vorlesung Mo 10:15 - 11:45, Physik PH 2501 (HS1)• Tutorubungen (vorauss.):

2x Mi 08:30 - 10:00 Uhr1x Mi 10:00 - 11:30 Uhr1x Mi 12:15 - 13:45 Uhr2x Fr 08:30 - 10:00 Uhr2x Fr 14:00 - 15:45 Uhr

Informationen• Webseite: www5.in.tum.de→ Teaching→ Einfuhrung in die

wissenschaftliche Programmierung• FRAGEN!




Organisatorisches II

Personen• Andreas Paul ([email protected]): Ubungsleitung• Benjamin Peherstorfer ([email protected]): Ubungsleitung• Christian, Johannes und Martin: Tutoren• Tobias Neckel ([email protected]): Vorlesung

Notige Software• Python 2.5.x, 2.6.x oder 2.7.x (idle, ipython);

ACHTUNG: Python 3.x nicht abwartskompatibel!• Bibliotheken: numpy, scipy, matplotlib• Beliebiger Editor (emacs, vi, nedit, ...)• Im CIP-Pool ist all das vorhanden




Organisatorisches IIIToDo• TUMOnline: Anmelden fur Veranstaltung• TUMOnline: Anmelden fur Tutorubung• TUMOnline: Anmelden fur Klausur

Bonus• Moeglichkeit 1 - 4 Punkte in der Klausur vorab zu erarbeiten• Voraussetzung: Anwesenheit in 9 unterschiedlichen Ubungen• Besprechung/Vorstellung einer (Teil-)Aufgabe in der Ubung

• Meldung direkt in Ubung (auf freiwilliger Basis)• Es werden 0 - 4 Punkte vergeben

• Im Fall von Krankheit• Falls Anwesenheit in 9 Ubungen nicht moglich, dann

arztliches Attest notig• Falls Vorstellung der Aufgabe zeitlich nicht mehr moglich,

dann schriftliche Bearbeitung einer Aufgabe am Ende




Typen von Programmiersprachen

Deklarative Programmierung• Funktional (Lisp, Scheme, make, ...)• Logisch (Prolog)• Sonstige (Sql, ...)

Imperative Programmierung• Prozedural (C, Fortran, (Python), Pascal, ...)• Modular (Bash, ...)• Objektorientiert (C++, Java, Python, ...)




Skripting vs. Compilieren

Skriptsprachen• Bash, Python, Perl, postscript, matlab/octave, ...• Interaktiver Modus• Teilweise wenig Optimierungen• Einfach zu verwenden• Keine Binardateien (daher fur kommerzielle Software meist

ungeeignet)Kompilierte Sprachen• C, C++, Fortran, ...• Effizient fur sehr rechenaufwandige Aufgaben• Sourcecode muss in Binarcode ubersetzt werden

Sonstige• Java




Warum Python? Warum Python in der Physik?

• Viele Physiker verwenden es (MPI etc.).• Pre- und Postprocessing sehr elegant!• Scripting!




Worum geht es in Python?

Python ist ...• High-level (Wirklich high-level!)• Komplett objektorientiert• Interpretiert• Skalierend• Erweiterbar• Portierbar (Windows, Linux, Mac OS X, Amiga, HPC, Cluster,

Web Server, Palm/Handhelds, Mobiltelefone, ...)• Vielseitig• Einfach zu lernen, zu verstehen, zu verwenden und zu warten• Kostenlos, open-source




Wenn ich nicht weiter weiß?

Literatur• RRZN Python-Skript beim LRZ (4,50 Euro)• David M. Beasley: Python - Essential Reference• Hans Petter Langtangen: A Primer on Scientific Programming

with Python• Vorlesungsfolien im Web• http://docs.python.org• www, google

Interaktive Hilfe• help() fur interactive Hilfe• help(’modules’) Liste der verfugbaren Module• help(object|’command’) Hilfe zu Objekt oder Befehl




Worum geht es in Python?

Kann verwendet werden fur:• OS Skripting (In vielen Fallen besser als Bash)• Internet Programmierung• Wissenschaftliches Rechnen (selbst komplexe Simulationen)• Parallelisierung• GUI (TKinter)• Visualisierung• Rapid Prototyping• ... und vieles mehr




Worum geht es in Python?Ursprunge• Erfunden 1990 als Lehrsprache• Hauptziel einfach lesbarer code• Benannt nach Monty Python

Weitere Features (nicht in C++!)• Dynamische Typisierung• Automatische garbage collection• Verschiedene Programmierparadigmen (Prozedural, OO,

funktional, Aspekt-orientiert, ...)Generelle Programmierregeln• Kommentare !• Problem⇒ Algorithmus⇒ Programm• Modulare Programmierung• Generisch wo moglich, spezifisch wo notig




Uberblick

Teil I: Erste SchritteTeil II: DatentypenTeil III: KontrollstrukturenTeil IV: Funktionen und ModuleTeil V: IO und Datentypen

Teil VI: Objektorientierte ProgrammierungTeil VII: Regulare AusdruckeTeil VIII: ExceptionsTeil IX: GrafikTeil X: PartikelsystemeTeil XI: Datenstrukturen

Teil XII: Wissenschaftliches RechnenTeil XIII: WarmeleitungTeil XIV: Losung von Linearen GleichungssystemenTeil XV: Mehrgitterverfahren




Python ausfuhren: python. . .

• Python interpreter starten (python oder ipython)• Python-prompt: >>> Befehl eintippen und mit <Enter>

bestatigen

$ python

Python 2.6.5 (r265 :79063 , Oct 1 2012, 22:04:36)

[GCC 4.4.3] on linux2

Type "help", "copyright", "credits" or "license" for

more information.

>>> help

Type help() for interactive help , or help(object) for

help about object.

>>> help()

Welcome to Python 2.6! This is the online help

utility. [...]




“Hello World!”Python – lesbarer Code• In Python

>>> print "Hello World!"

Hello World!

• Andere Beispiele: Java

public class Hello

public static void main(String argv [])

System.out.println("Hello World!");

$ javac Hello.java

$ java Hello

Hello World!




Variablen, Zuweisungen und Ausdrucke

>>> 3 + 4**2/8 + 1

>>> (6*3.5)*2 - 5

>>> 10e-4

>>> a = 17

>>> a

>>> a - 1.5

>>> a = a - 1

>>> b = "Hello World!"

>>> print b

>>> print a, b, (4 + 5**0.5)

• Zuweisung: Variablenname = Ausdruck• Nicht zu verwechseln mit mathematischem =• Ausgabe des Ergebnisses nur interaktiv• print gibt eine String-Reprasentation von Objekten aus• Leerzeichen im Ausdruck egal (Aber nicht VOR dem Ausdruck)




Variablennamen

• a-z, A-Z, 0-9,• Erstes Zeichen nicht 0-9• ist erlaubt, hat aber spezielle Bedeutung• Variablennamen sollten aussagekraftig sein• Bei mathematischen Ausdrucken entsprechend den

mathematischen Variablen• Ansonsten deskriptiv• Einige Worte sind reserviert:

and, as, assert, break, class, continue, def, del, elif, else,

except, False, finally, for, from, global, if, import, in , is,

lambda, None, nonlocal, not, or, pass, raise, return, True,

try, with, while, yield

• Variablen in Python sind dynamisch typisiert




3.5 Wege zu Python

1. Python interpreter starten mittels python

2. IPython nutzen: ipython3. Python auf einer Datei ausfuhren

3.5. Python uber ein ausfuhrbares Skript nutzen




2. Python ausfuhren: ipythonVorteile• Erweiterte interaktive Shell• Zusatzliche Shell Syntax

• “magische” Funktionen, beginnend mit “%”, z.B. %psearchum Funktionen im Namensraum zu suchen

• Zugriff auf Bash-Befehle %cd

• Python-Skripte ausfuhren mit run filename.py

• Logging, Makros, . . .• Pretty printing, umschalten %Pprint

• Syntax highlighting• Tab-Vervollstandigung• History

• Vorherige Befehlsblocke wiederholen• Automatische Einruckung• Mitgeliefert mit SciPy




3. Python ausfuhren: Datei ubergeben

• Datei square_me.py:

print 6**2

• Ausfuhren mit

$ python square_me.py

36




3.5 Python ausfuhren: ausfuhrbares Skript

• Gleicher Programmcode wie zuvor• Dem Betriebssystem mitteilen, dass python verwendet werden

soll:

#!/usr/bin/python

print 6**2

• Datei ausfuhrbar machen und ausfuhren

$ chmod u+x square_me.py

$ ./ square_me.py

36




Teil II

Datentypen




Numerische Typen - ganze Zahlen (int)

• Rechnen, wie wir es gewohnt sind (12*34+567)• In Python ist der Zahlenbereich beliebig groß

>>> 1234**56

12991190255487145194103208439623513775465782010127

39238437901270462425943305509464892567848536247290

20106139515647384910944921186523865849056275359066

262352911682504769929216L

>>> type (1234)

<type ’int’>

>>> type (1234**56)

<type ’long’>




Numerische Typen - float und bool

• bool – wahr und falsch

>>> print True , False

>>> True == True

>>> True != True

• float – Gleitkommazahl (double precision)f = s ·m · 2e, 1+52+11 Bit

>>> 1.0 - 49.0*(1.0/49.0)

1.1102230246251565e-16

>>> 1 - 49*(1/49)

1

>>> 1 - 49*(1/49.0)

1.1102230246251565e-16

• None – special type (“undefined”)




Numerische Typen - Komplexe Zahlen

>>> 1+4j

(1+4j)

>>> complex (2,3)

(2+3j)

>>> 2*(c**3 - 3j)

( -92+12j)

>>> d.real

-92.0

>>> d.imag

12.0

• Imaginarteil hat Zusatz j

• Basisoperationen normal anwendbar• Fur komplexere Operationen (z.B. sin) spezielle Bibliotheken




Allgemeine Hinweise• Alles in Python ist ein Objekt! (Attribute und Methoden)• Zugriff auf Objektmethoden: <Objektname>.<Methodenname>()• Zugriff auf Objektattribute: <Objektname>.<Attributname>

c = complex (2,3)

print c.real , c.imag

• Dynamische Typisierung: Typ von Variablen durch Zuweisung fix• Automatische Konvertierung wo notig und moglich

a = 1234

type(a)

b = 100

type(b)

a = 1234**100

type(a)

b = b*1.0 # entspricht b = float(b)*1.0

type(b)




Numerische Typen - Arithmetische Operationen

• Klassische Operationen

x + y # Addition

x - y # Subtraktion

x * y # Multiplikation

x / y # Division

x // y # Truncated division (integer division)

x ** y # Exponentiation

x % y # Modulo

x -= 2; x *= 4; ...

• Achtung: Division verschieden fur int und float (Python < 3.0)

7/4

7.0/4




Numerische Typen - Vergleichsoperatoren

• a == b gibt True bei Gleichheit (= dient Zuweisung)

>>> 1 == 2

False

>>> 1 == 1.0

True

>>> 1 == "1"

False

• weitere Operatoren: !=, <, >, <=, >=

• Verknupfung mit and und or, Negation mit not• Shifting: <<, >>

• Bitweise Operatoren: &, |, ^, ~

>>> not(1 == 2 or 3 + 1 == 2 << 1)

False




Mathematische Funktionen: Modul math

• Ein Modul ist eine Sammlung von Funktionalitaten• Einbinden mit import• math enthalt mathematische Standardfunktionen• sin, cos, sqrt, exp, log, pow, ...

• Zusatzlich Konstanten pi und e

• Funktioniert nur mit integer und float• Fur komplexe Zahlen: cmath

>>> import math

>>> math.sqrt (2)

1.4142135623730951

>>> from math import sin

>>> sin (1)

0.8414709848078965




Sequenzen - Strings• Sequenz von Zeichen (character)• einzelne oder doppelte Anfuhrungszeichen• dreifache Anfuhrungszeichen fur mehrzeiligen String

"... there lived a hobbit."

"a ’hobbit ’"

’a "hobbit"’

’a \’hobbit\’’

"""In a hole in the ground

there "lived" a ’hobbit ’"""

• Verknupfung von Strings

s1 = "In a hole in the ground"

s2 = "there lived a hobbit"

s = s1 + ’ ’ + s2

print s




Sequenzen - Strings (2)

• Replikation

>>> "Hi! "*3

’Hi! Hi! Hi! ’

• Indizierung und Abschnitte

>>> s = "12345678"

>>> s[3], s[2:4], s[-2]

(’4’, ’34’, ’7’)

>>> s[6:], s[:3]

(’78’, ’123’)

>>> s[0: -1:2] # stride 2

’1357’

>>> s[2] = 4 # Fehler: Strings sind unveranderbar

>>> len(s)

8




Ausgewahlte String-Methoden

s = " Frodo and Sam and Bilbo"

s.islower ()

s.isupper ()

s.startswith("Frodo") # s.startswith (" Frodo", 2)

s.endswith("Bilbo")

s = s.strip()

s.upper ()

s = s.lower()

s = s.center(len(s)+4)# zentrieren

s.lstrip ()

s.rstrip(" ")

s = s.strip()

s.find("sam")

s.rfind("and")




Ausgabe verschiedener Datentypen

print

>>> a=3; b=4.5; c=’text’

>>> print a, b, c

3 4.5 text

Stringinterpolation• Stringinterpolation gibt hohe Flexibilitat• Platzhalter mit % in einem String• Nach dem String folgen ein % und ein Tupel mit Werten.

>>> print "int: %d, float: %f und string: %s" \

% (4, 3.5, "bla")

int: 4, float: 3.500000 und string: bla




Konvertierungszeichen zur Stringinterpolation

• %d Dezimal Integer• %f Gleitpunkt (float)• %e,%E Gleitpunkt wissenschaftlich (m.ddde+xx)• %g,%G kompakteste Gleitpunktdarstellung• %s String• %c einzelnes Zeichen• %o Oktaldarstellung eines Integers• %x, %X Hexadezimaldarstellung eines Integers• %% Das Zeichen %

>>> a = 7.357

>>> print "%f, %.2f, %8.1f" % (a, a, a)

7.357000 , 7.36, 7.4




Sequenzen - Tupel

• Sequenzen beliebiger Objekte• Unveranderlich• Indizierung und Ausschnitte wie bei Strings

>>> t = (1, "zwei", (3,4,5))

>>> t[0]

1

>>> t[-1]

(3,4,5)

>>> t[1] = 2 # error

>>> len(t)

3

>>> t + (6,7)

(1, ’zwei’, (3, 4, 5), 6, 7)

• Minimun und Maximum mit min(t) und max(t)




Sequenzen - Listen• Sequenzen beliebiger Objekte• Veranderlich• Indizierung und Ausschnitte wie bei Strings

>>> l = [1, "zwei", (3,4,5)]

>>> l[0]

1

>>> l[1] = 2; l # funktioniert!

[1, 2, (3,4,5)]

>>> l + [6,7]

[1 , 2, (3, 4, 5), 6, 7]

• Operationen auf Listen

l = [0,1,2,3,4,5,6,7,8,9]

l[2:4] = [10, 11]

l[1:7:2] = [-1, -2, -3]

del l[::2]




Sequenzen - Listen (2)• Listen sind Objekte mit Methoden:

>>> l = [0,1,2]

>>> l.append("x") # [0, 1, 2, ’x’]

>>> l.extend ([5 ,6]) # [0, 1, 2, ’x’, 5, 6]

>>> l.insert(2, "x") # [0, 1, ’x’, 2, ’x’, 5, 6]

>>> l.count("x") # 2

>>> l.sort() # [0, 1, 2, 5, 6, ’x’, ’x’]

>>> l.reverse () # [’x’, ’x’, 6, 5, 2, 1, 0]

>>> l.remove("x") # [’x’, 6, 5, 2, 1, 0]

>>> l.pop() # [’x’, 6, 5, 2, 1]

0

• Integerlisten erzeugen: range Funktion

>>> range (5)

[0, 1, 2, 3, 4]

>>> range(-10, 5, 3)

[-10, -7, -4, -1, 2]T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Teil III

Kontrollstrukturen




Struktogramme

• Graphische Darstellung von Programmen• Anweisungen

Anweisung 1Anweisung 2

...Anweisung n

• BedingungenBedingung

wahr falschAnweisungs-

block 1Anweisungs-

block 2

• SchleifenSchleifenbedingung

Anweisungsblock




Struktogramme - Beispiel

• Euklidscher Algorithmus zur Berechnung des ggT:Wenn CD aber AB nicht misst, und man nimmt bei AB, CDabwechselnd immer das kleinere vom großeren weg, dann muss(schließlich) eine Zahl ubrig bleiben, die die vorangehende misst.

[Euklid, Die Elemente, Buch VII, §2]

Solange a>0 und b>0

a = a - b

a>b?ja nein

b = b - a

Ausgabe a

b==0?ja nein

Ausgabe b




KontrollstrukturenGenerelle Bemerkungen• Blocke (siehe Struktogramme) mussen im Code markiert werden• In vielen Sprachen ublich: geschweifte Klammern• Python verwendet stattdessen Einruckungen• Einrucktiefe entscheidet uber Blockzugehorigkeit• Doppelpunkt leitet Block ein

Fallunterscheidung: if

if x < y:

print "x is smaller"

print x

elif x == y:

print "both are equal"

else:

print "y is smaller"

print y




Nochmal Typen und Objekte• Jedes Objekt hat Identitat (id), Typ (type) und Wert.

>>> b = 42 # Wert: 42

>>> type(b)

<type ’int’>

>>> id(b)

158788940

>>> type(type(b))

<type ’type’>

• Vergleich zweier Objekte:

if a is b:

print ’a und b sind das identische Objekt ’

if a == b:

print ’a und b haben denselben Wert’

if type(a) is type(b):

print ’a und b haben denselben Typ’




Nochmal Typen und Objekte (2)

• Wenn zwei Objekte identisch sind (is bzw. id(a) == id(b)),haben sie auch den gleichen Wert und Typ.

• Gleicher Typ sagt nichts uber Identitat und Wert aus.• Gleicher Wert sagt nichts uber Identitat aus.

>>> b = 42

>>> l = [1,2,3]

>>> type(b) is list

False

>>> type(l) is list

True

>>> filehandle = open("test.txt", "w")

>>> type(filehandle)

<type ’file’>




Kontrollstrukturen - Zahlschleife

for iteriert uber die Elemente einer Sequenz.

for i in [1,2,3]:

print "Zahl: ", i

for i in "Hallo":

print "Buchstabe: ", i

for i in (1, ’a’, [4,3,2], "Hallo"):

print "Element: ", i




Kontrollstrukturen - Schleifen mit Bedingungen• continue startet den nachsten Schleifendurchlauf (for und while).• break bricht eine Schleife ab (for und while).

a = ("let", "us", "find", "Bilbo", "again")

for word in a:

if word == "Bilbo":

print "found Bilbo"

break

print word

• Eine while Schleife iteriert, solange ihre Bedingung True ist.

a = 2048; a_orig = a; exp = 0

while a > 1:

a /= 2 # a = a / 2

exp+=1 # exp = exp + 1

print a_orig , "=", 2, "^", exp




Euklidscher Algorithmus mit Python

while a>0 and b>0:

if a>b:

a = a-b

else:

b = b-a

if (b==0):

print a

else:

print b




Beispiel: Numerische Losung von ODEs• Logistische Differentialgleichung: p(t) = ap(t)− bp(t)2 (Modell

von Verhulst)• Anfangswert p(0) = p0• Losung p(t) = a·p0

b·p0+(a−b·p0)·e−at

• Z.B. a = 1, b = 1/100

p(t)

t0 1 2 3 4 5 6 7 8 90

20

40

60

80

100




Analytische LosungMittels einer Schleife konnen wir uns eine Wertetabelle druckenlassen:

from math import exp

a = 1. # Parameter der DGL

b = 0.01

p0 = 10. # Anfangswert

tend = 5. # Intervall [0, tend]

N = 10 # Teilintervalle

for i in range(0, N+1):

t = (tend*i)/N

p = a*p0 \

/(b*p0 + (a-b*p0)*exp(-a*t))

print ’%6.1f %6.1f’ % (t, p)

0.0 10.0

0.5 15.5

1.0 23.2

1.5 33.2

2.0 45.1

2.5 57.5

3.0 69.1

3.5 78.6

4.0 85.8

4.5 90.9

5.0 94.3

Was passiert, wenn wir in tend = 5. den Punkt weglassen?




Numerische Losung: Explizites Euler-Verfahren• Verfahren zur numerischen Losung einer DGL• DGL: p(t) = f (t ,p(t))

• Ersetze Ableitung durch Vorwartsdifferenz p(t) = p(t+δt)−p(t)δt

• ⇒ p(t + δt)=p(t) + δt · f (t ,p(t))

• δt positiv und betragsmaßig klein• Berechnungsvorschrift, um von p(t) einen Schritt δt in die

Zukunft zu gehen




Das Euler-Verfahren fur die logistische Differentialgleichung, vierverschiedene Zeitschrittweiten δt :

0

10

20

30

40

50

60

0 0.5 1 1.5 2 2.5 3

p(t)

t

Loesung der logistischen DGL, a=1, b=0.01, p0=10Euler-Verfahren, dt=1.000Euler-Verfahren, dt=0.500Euler-Verfahren, dt=0.250Euler-Verfahren, dt=0.125




Das Programm fur das Euler-Verfahren sieht ganz ahnlich aus wiedas fur die Wertetabelle von vorhin:


b = 0.01


tend = 3. # Interv. [0, tend]

N = 6 # Teilintervalle

dt = tend/N # Zeitschrittweite

p = p0

t = 0.

for i in range(0, N):

t = t + dt

p = p + dt*(a*p - b*p**2)

print ’%6.1f %6.1f’ % (t, p)

0.5 14.5

1.0 20.7

1.5 28.9

2.0 39.2

2.5 51.1

3.0 63.6




Zusammenhang zwischen Zeitschrittweite und Genauigkeit• Programm umbauen, um in Abhangigkeit von der Arbeit

(Teilintervalle N) die Genauigkeit (Fehler) anzugeben• Vergleich nur am Ende des Intervalls• from math import exp zur Berechnung der exakten Losung• pend = p0*a/(p0*b + (a-b*p0)*exp(-a*tend))

• print ’%6d %7.3f %7.3f’% (N, p, pend-p)

• Fur verschiedene N = 6 · 2j , j = 0...9:

for j in range (0 ,10):

N = 6*2**j # Anzahl Teilintervalle




Erweitertes Programm

>from math import exp


b = 0.01


tend = 3. # Interv. [0, tend]

>pend = p0*a/(p0*b + (a-b*p0)*exp(-a*tend))

>for j in range (0 ,10):

>N = 6*2**j # Anzahl Teilintervalle

dt = tend/N # Zeitschrittweite

p = p0

t = 0.

for i in range(0, N):

t = t + dt

p = p + dt*(a*p - b*p**2)

>print ’%6d %7.3f %7.3f’ % (N, p, pend -p)




Ergebnisse EulerN Naherung Fehler6 63.590 5.467

12 66.529 2.52824 67.847 1.20948 68.466 0.59196 68.765 0.292

192 68.912 0.145384 68.984 0.072768 69.021 0.036

1536 69.039 0.0183072 69.048 0.009

• Halbierung von δt (Verdoppelung von N) halbiert den Fehler• Vgl. Folie 51




Numerische Losung: Verfahren von Heun• Konvergenz des Euler-Verfahrens schlecht• Idee zur Verbesserung: Nicht nur Steigung am aktuellen Punkt• Mittelwert aus Steigung am aktuellen Punkt und Steigung am

nachsten Punkt (bzw. dessen Naherung):

p(t + δt) .= p(t) + δt

f (t ,p(t)) + f (t + δt ,p(t) + δt · f (t ,p(t)))

2

• bisherige Programmzeile: p = p + dt*(a*p - b*p**2)

• wird ersetzt durch;

f1 = a*p - b*p**2

ptmp = p + dt*f1

f2 = a*ptmp - b*ptmp **2

p = p + dt*(f1 + f2)/2




Das Heun-Verfahren fur die logistische Differentialgleichung, dreiverschiedene Zeitschrittweiten δt :

0

10

20

30

40

50

60

0 0.5 1 1.5 2 2.5 3

p(t)

t

Loesung der logistischen DGL, a=1, b=0.01, p0=10Verfahren von Heun, dt=1.000Verfahren von Heun, dt=0.500Verfahren von Heun, dt=0.250




Ergebnisse HeunN Naherung Fehler6 68.140 0.916

12 68.807 0.25024 68.992 0.06548 69.040 0.01796 69.053 0.004

192 69.056 0.001

• Halbierung von δt (Verdoppelung von N) viertelt den Fehler• Euler und Heun sind Einschrittverfahren• Neben dem Fehler sind weitere Eigenschaften wichtig (z.B.

Energieerhaltung)• Komplexere Einschrittverfahren, z.B. Runge-Kutta• Daruber hinaus Mehrschrittverfahren




Ausblick

Was fehlt?• Kommentare• An einigen Stellen wurde Code dupliziert• ⇒ schlechte Wartbarkeit• Irgendwann wird der Code lang und unubersichtlich

Losung?• Codestucke zusammenfassen und auslagern in Funktionen• Ubergabe von Parametern an die Funktion• Auswertung der Ruckgabewerte der Funktion




Teil IV

Funktionen und Module




Wozu Funktionen?Wo sie schon verwendet wurden:• Mathematische Funktionen: sqrt(s), sin(s), exp(x), ...

• Methoden von Sequenzen: s.islower(), l.append(x), ...

• Erzeugung von Listen: range(x)• Typ und Identitat: type(x), id(x)

• Alle bisherigen haben keinen oder einen Ubergabewert undeinen Ruckgabewert

Vorteile• Zerlegung eines großen Problems in viele kleine• Zusammenfassung von Anweisungsblocken• Klare Abgrenzung von Funktionalitaten• Vermeidung von Code-Duplikation• Hohere Code-Lesbarkeit (Wenn man statt sin(x) jedesmal den

kompletten Code ...)• ...




Syntax von Funktionen• Funktionen werden definiert mit dem Schlusselwort def• Nach Funktionsname und Parameterliste kommt ein

Doppelpunkt.• Der darauf folgende Block (Einruckung!) wird bei Aufruf der

Funktion ausgefuhrt.• Ausfuhrung bis zum Ende des Blocks oder bis zu einen return

• Ruckgabewert kann beliebiger Typ sein (Zahl, Liste, String, ...)• Ohne return oder ohne Wert wird None zuruckgegeben

>>> def hi():

>>> print "Hello World!"

>>>

>>> hi()

Hello World

>>> a = hi()

Hello World

>>> type(a)

<type ’NoneType ’>T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Funktionen sind Objekte

>>> def factorial(n):

>>> fac = 1

>>> for i in range(2, n+1):

>>> fac = fac*i

>>> return fac

>>>

>>> factorial (3)

6

>>> f = factorial

>>> f(4)

24

>>> type(f)

<type ’function ’>




Funktionsparameter

• Parameter stehen in runden Klammern hinter demFunktionsnamen.

• Formalparameter: Bei der Definition der Funktion• Aktualparameter: Folge von Ausdrucken beim Aufruf der

Funktion• Beliebig viele Parameter moglich, normalerweise gleich viele

Formal- und Aktualparameter

>>> def potenziere(basis , exponent ):

>>> return basis ** exponent

>>>

>>> potenziere (2,10)

1024




Standardwerte fur Parameter• Formalparameter konnen mit Standardwerten belegt werden.• Zuerst alle Formalparameter ohne Standardwerte, dann die mit

def potenziere(basis=2, exponent ): # Falsch

def potenziere(basis , exponent =10):

def potenziere(basis=2, exponent =10):

>>> potenziere (2 ,10) # (2 ,10)

>>> potenziere (2) # (2,2)

>>> potenziere(exponent =7) # (2,7)

• Fur diese Funktion sind Standardwerte also relativ nutzlos.• Bei vielen Parametern oder bestimmten Situationen kann es

aber hilfreich sein.




Ruckgabewerte• Bei mehreren Ruckgabewerten wird ein Tupel zuruckgegeben.

>>> def return_two ():

>>> return "x", "y"

>>>

>>> h = return_two (); # h = (’x’, ’y’)

>>> (a,b) = return_two () # a = ’x’, b = ’y’

• Eine Liste oder explizites Tupel ist ein einzelner Wert.

>>> def primfaktoren(x):

>>> l=[]; n=2

>>> while n <= x:

>>> while x%n==0:

>>> x /= n

>>> l.append(n)

>>> n += 1

>>> return l




Gultigkeitsbereich

Funktionen• Erinnerung: Python ist eine interpretierte Sprache!• Funktionen mussen definiert sein, bevor sie aufgerufen werden

konnen.Variablen• Unterscheidung zwischen lokalen und globalen Variablen• Eine Variable, die in einer Funktion definiert (an einen Wert

gebunden) wird, ist lokal und insbesondere außerhalb nichtsichtbar.

• Grundregel: GLOBALE VARIABLEN VERMEIDEN!!!• Globale Variablen konnen gelesen werden.• Um globale Variablen in einer Funktion zu schreiben, wird das

Schusselwort global verwendet.




Beispiele

>>> def zaehle ():

>>> global anzahl

>>> anzahl += 1

>>>

>>> anzahl = 0

>>> zaehle (); zaehle (); zaehle ()

>>> anzahl

3

>>> def no_effect ():

>>> xyz = 5**2

>>>

>>> no_effect ()

>>> xyz

NameError: name ’xyz’ is not defined




Docstrings• Ein Grund fur Funktionen: Funktionalitat kapseln• Jemand, der die Funktion nicht geschrieben hat, soll sie

verwenden konnen.⇒ Dokumentation notig!• Ein String (ein- oder mehrzeilig) nach dem Header ist fur

interaktive Dokumentation vorgesehen.• Python ignoriert den String (Ausdruck ohne Zuweisung) bei der

Funktionsausfuhrung.• Das Hilfesystem (und der Code-Leser) kann den String

verwenden.• Nicht nur fur Funktionen, sondern auch fur Klassen, Module, ...

>>> def blubb ():

>>> "Unsinnige Funktion , die ’blubb ’ ausgibt"

>>> print "blubb"

>>> blubb()

blubb

>>> help(blubb)




Call by Object Reference

• Bei Funktionsaufruf: Referenz auf Parameterobjekt (”Adresse“ ,id) wird ubergeben

• Referenz (id) wird kopiert (Vorteil: wenig Daten)

Object

3

Ref.

X

in fct.

X

Ref.

outside




Call by Object Reference (2)• ACHTUNG: Zuweisung andert Referenzobjekt in der Funktion!

>>> def add_one(x):

>>> x = x + 1

>>> x = 3

>>> add_one(x)

>>> x #still x==3

Object

3X

Ref.

outside outside

Ref.

X

in fct. in fct.

Object

4




Call by Object Reference (3)

• Nach der Zuweisung zeigt die lokale Variable auf ein neuesObjekt.

• Mit Methoden kann dieser Effekt verhindert und das bisherigeObjekt verandert werden.

>>> def append1(l):

>>> l = l + [4]

>>>

>>> def append2(l):

>>> l.append (4)

>>>

>>> l = [1,2,3]

>>> append1(l); l

[1,2,3]

>>> append2(l); l

[1,2,3,4]




Anonyme Funktionen

• Normale Funktionen mussen separat definiert werden.• Sie bekommen einen Namen, mit dem auf sie zugegriffen wird.

>>> def inc(i):

>>> return i+1

>>> inc (3)

4

• Namenlose/Anonyme Funktionen erhalten keinen Namen.• Sie konnen an beliebigen Stellen in den Code eingebaut werden.

>>> inc = lambda(i): i+1

>>> inc (3)

4




map• map wendet eine Funktion auf alle Elemente einer Liste an• Die Funktion kann naturlich auch anonym sein.

>>> L = range (10)

>>> map(lambda x: x*x+1, L)

[1, 2, 5, 10, 17, 26, 37, 50, 65, 82]

• In diesem Fall ware das sogar noch kurzer mit listcomprehensions (nachste Folie) gegangen.

filter• filter wendet ebenfalls eine Funktion auf alle Elemente einer

Liste an.• Diejenigen Elemente, fur die die Funktion True zuruckgibt,

werden zuruckgegeben.

>>> filter(lambda x: x%2==0, L)

[0, 2, 4, 6, 8]




List Comprehension

• Einfache Listen lassen sich z.B. mit range erzeugen.• Fur komplexere Listen gibt es list comprehension.

>>> b = [i**2 for i in range (0 ,10) if i%2==0]

>>> b

[0, 4, 16, 36, 64]

• list comprehensions konnen geschachtelte Schleifen nachbilden:

>>> c = [(i,j) for i in range (1,5) if i%2==0

for j in range (1,5) if j%2!=0]

>>> c

[(2, 1), (2, 3), (4, 1), (4, 3)]




List Comprehension (2)• Generelle Syntax:

[expression for item1 in iterable1 if condition1

for item2 in iterable2 if condition2

...

for itemN in iterableN if conditionN]

• Falls keine Bedingung angegeben: automatisch True

• Syntax entspricht:

s = []

for item1 in iterable1:

if condition1:

for item2 in iterable2:

if condition2:

....

for itemN in iterableN:

if conditionN: s.append(expression)




Beispiel: Numerische Integration (Quadratur)Berechnung des Integrals

F1(f ,a,b) :=

∫ b

af (x) dx

Verschiedene numerische Verfahren; hier jetzt:Trapezregel• Bestimmung des Funktionswertes an beiden Randern• Multiplikation des Mittelwerts mit der Intervallbreite

• F1 ≈ T := (b − a) · f (a)+f (b)2

• Polynom 1. Ordnung (Flache entspricht Trapez)Quadraturfehler

|T − F1| ≤112· sup

x∈[a,b]|f ′′(x)| · (b − a)3




Summenregeln• Bei großen Intervallen wird auch der Fehler sehr groß.• ⇒ Unterteilung in n kleine Intervalle• Jedes Teilintervall hat Lange h = (b − a)/n

Trapezsumme

TS := h ·[

f (a)

2+

n−1∑

i=1

f (a + ih) +f (b)

6

]




Implementierung der Trapezsumme

def trapezsumme(f, a, b, n):

h = (b-a)/n

sum = (f(a) + f(b))/2

for i in range(1,n):

sum = sum + f(a + h*i)

return sum*h

Funktion zum Test• Funktion, die analytisch leicht zu integrieren ist• ⇒ Fehler lasst sich leicht berechnen

import math

def f(x):

return math.sin(math.pi*x)

def f_int(x):

return -math.cos(math.pi*x)/math.pi




Berechnung mit unterschiedlicher Genauigkeit: Fehler empirisch

• n = 21...26

a_bsp = 0.

b_bsp = 1.

exakt = f_int(b_bsp) - f_int(a_bsp)

for i in range (1,7):

n = 2**i

t = trapezsumme(f, a_bsp , b_bsp , n)

print ’%4d % -12.6g % -12.6g’ % (n, t, exakt -t)

n Naherung Fehler2 0.5 0.136624 0.603553 0.03306648 0.628417 0.00820234

16 0.634573 0.0020466232 0.636108 0.00051140964 0.636492 0.000127837




Fehler der Trapezsumme: Fehler analytisch• In jedem Teilintervall ist der Fehler

≤ 112· sup

x∈[a,b]|f ′′(x)| · h3

• Bei n = (b − a)/h Intervallen ist damit der Gesamtfehler

|TS − F1| ≤112· sup

x∈[a,b]|f ′′(x)| · (b − a) · h2

• Verdopplung des Rechenaufwands (Halbierung von h) reduziertden maximalen Fehler um den Faktor 4




Module• Funktionalitat zusammenfassen in einer separaten .py Datei• Importieren Bibliotheksmodul (Schon verwendet: math, time, ...)• Beispiel (tools.py):

""" This module provides some helper tools.

Try it."""

counter = 42

def readfile(fname):

"Read text file. Returns list of lines"

fd = open(fname , ’r’)

data = fd.readlines ()

fd.close ()

return data

def do_nothing ():

"Do really nothing"

pass




Importieren• Modul importieren

import tools

tools.do_nothing ()

print tools.counter

• Modul unter neuem Namen importieren

import tools as t

t.do_nothing ()

• Einzelne Funktionen direkt importieren

from tools import do_nothing , readfile

from tools import counter as cntr

do_nothing ()

print cntr




Import beeinflussen• Alle Symbole in den Namensraum importieren

from tools import *

do_nothing ()

print counter

• Module konnen bestimmen, welche Symbole mitfrom module import * importiert werden:

# module tools.py

__all__ = [’readfile ’, ’counter ’]

Die Funktion do_nothing() ist nach import * nicht bekannt!Funktionalitat abfragen•

dir(tools)

dir(math)




Hilfe zu Modulen und Funktionen• Docstrings abfragen mit help

import tools

help(tools)

help(tools.do_nothing)

Ausfuhrbares Programm als Modul?• tools.py soll sowohl Bibliothek sein, als auch selbst ausfuhrbar

# tools.py

...

if __name__ == ’__main__ ’:

print "tools.py executed"

else:

print "tools.py imported as module"




Module debuggen• Wenn ein Modul geandert wird, werden Anderungen nur

wirksam, wenn Python neu gestartet wird, oder reload(Python<3.0) verwendet wird.

reload(tools)

Pfad fur Module• Python sucht im Suchpfad und im aktuellen Verzeichnis• Suchpfad kann erweitert werden

import sys

sys.path.append("folder/to/module")

import ...




Pakete

• Module konnen gruppiert werden• Hierbei ist die Verzeichnisstruktur wichtig

tools/

__init__.py # Inhalt fur "import tools"

files.py # fur "import tools.files"

graphics.py # fur "import tools.graphics"

stringtools/

__init__.py # fur "import tools.stringtools"

... weitere Verschachtelung moglich

• Wenn from tools import * Untermodule enthalten soll, musstools/__init__.py folgende Zeile enthalten:

__all__ = ["files", "graphics"]




Teil V

IO und weitere Datentypen




Datei Ein- und Ausgabe

Dateien offnen• Textuelles Datei-Objekt erstellen

datei = open("testfile.txt") # read

datei = open("testfile.txt", ’r’) # read

datei = open("testfile.txt", ’w’) # write

datei = open("testfile.txt", ’a’) # append

• Binares Datei-Objekt erstellen

datei = open("testfile.txt", ’rb’) # read

datei = open("testfile.txt", ’wb’) # write

datei = open("testfile.txt", ’ab’) # append

• open Hat noch weitere Optionen (Codierung, Behandlung vonnewlines, . . . )




Methoden von Datei-Objekten• Lesen

datei.read() # komplett einlesen

datei.read(n) # n Byte einlesen

datei.readline () # eine Zeile lesen

datei.readlines () # Liste von Zeilen

• Schreiben

datei.write("Neuer Text\n")

datei.writelines (["Erste Zeile\n", "Zweite Zeile\n"])

• Nicht vergessen: newlines (\n)• Datei schließen

datei.close ()




Textdateien zeilenweise bearbeiten• for Schleife uber die Zeilen der Datei

datei = open("fulltext.txt")

for line in datei: # alt.: in datei.readlines ()

print line

# oder:

while True:

line = datei.readline ()

if not line:

break

print line

Weitere Methoden von Dateiobjekten• datei.tell() – aktuelle Position ausgeben• datei.seek(pos) – Zur gegebenen Position gehen• datei.flush() – Ausgabepuffer leeren (Pufferung wegen

Effizienz)




Standard Input, Output, und Error

• Modul sys stellt drei Standard-Dateiobjekte bereit• sys.stdin – nur lesen• sys.stdout, sys.stderr – nur schreiben

import sys

sys.stdout.write("Text eingeben: ") # = print "..."

line = sys.stdin.readline ()

# oder: line = raw_input ("Text eingeben: ")

if error:

sys.stderr.write("Error!\n")

• write Hangt nicht automatisch einen Zeilenumbruch an• raw_input([text]) entfernt Zeilenumbruche• Ein- und Ausgabe sind gepuffert




Kommandozeilenparametersys.argv• sys.argv ist eine Liste mit Kommandozeilenparametern• Erstes Element (sys.argv[0]) ist der Programmname

import sys

print "Programmname: %s" % (sys.argv [0])

if len(sys.argv) < 2:

print "Keine Parameter angegeben"

sys.exit (1)

print "Parameters:"

print ", ".join(sys.argv [1:])

• Bei sehr wenigen Parametern ist diese Methode praktikabel• Schlecht bei vielen Parametern, oder wenn Parameter nur

optional sein sollen⇒ Modul argparse




nochmal Sequenzen: Dictionaries und Sets

Bisherige Container zur Speicherung großerer Datenmengen• Tupel: unveranderlich• Liste: veranderlich• Beide werden indiziert mit 0...n-1 (n Anzahl an Elementen)

Vorteile• Tupel: sehr effizient, da unveranderlich• Liste: hoher Overhead, aber immer noch effizient• Optimal, wenn eine Folge von Werten gespeichert werden soll

Beschrankungen• Keine Mengen darstellbar• Keine Bezeichner fur Objekte moglich

z.B. suche Buch mit bestimmter ISBN oder bestimmtem Titel




Dictionary• Bildet einen Schlussel auf einen Wert ab• Werte konnen beliebigen Typ haben• Schlussel nur unveranderliche Objekte (int, float, string, tuple,...)• Kann als primitive Datenbank verwendet werden• Erstellung mit geschweiften Klammern• Optional in den Klammern: Folge von key:value Paaren

>>> nummern =

>>> nummern["Tobias"] = "089 289 18632"

>>> nummern["Benjamin"] = "089 289 16830"

>>> nummern["Benjamin"]

"089 289 18630"




Methoden von Dictionaries

>>> d = 7:’A’, ’Muh’: ’Maeh’, (4, 5, 6):[6, 5, 4]

>>> d.keys()

[(4, 5, 6), ’Muh’, 7]

>>> d.values ()

[[6, 5, 4], ’Maeh’, ’A’]

>>> d.items()

[((4, 5, 6), [6, 5, 4]), (’Muh’, ’Maeh’), (7, ’A’)]

>>> del d[’Muh’]

>>> for (key , value) in d.items ():

>>> print key , value

>>> for i in d.items ():

>>> print i[0], i[1]

>>> for key in d.keys ():

>>> print key , d[key]




Fehlerhafte Zugriffe• Tupel und Listen: gesamter Indexbereich mit Werten belegt• Bei Dictionaries gibt es keinen Indexbereich• Was passiert, wenn ein unbekannter Schlussel gesucht wird?

>>> nummern[’Andreas ’]

KeyError: 4

>>> nummern.has_key(’Benjamin ’)

True

>>> ’Benjamin ’ in nummern

True

Standardwerte• Methode get liefert None bei nicht existierendem Schlussel• Optionaler zweiter Parameter fur get definiert Standartwert

>>> nummern.get(’Benjamin ’)

"089 289 18630"

>>> nummern.get(’Andreas ’, ’keine Nummer vorhanden ’)

’keine Nummer vorhanden ’




Set• Reprasentiert eine Menge• Jedes Element ist maximal ein Mal in der Menge• set ist veranderbar, Elemente mussen aber unveranderlich sein• Kann initialisiert werden mit beliebiger Sequenz

set("Hallo") # set([’a’, ’B’, ’l ’])

s=set([3,1,2,3,"Bla" ,2]) # set([1, 2, 3, ’Bla ’])

l = [2,8,7]

s.difference(l) # set([1, 3, ’Bla ’])

s.intersection(l) # set ([2])

s.union(l) # set([1, 2, 3, 7, 8, ’Bla ’])

s.symmetric_difference(l) # set([1, 3, 7, 8, ’Bla ’])

s.issubset ([3,"Bla"]) # False

s.issuperset ([3,"Bla"]) # True

• Weitere Methoden: add, remove, ...

• bestimmte Operatoren auf sets unterstutzt: -, ˆ, |, &




Interaktion mit dem BetriebssystemModul os• Schnittstelle zu Funktionalitaten des Betriebssystems• Variablen:

import os

os.environ # Dictionary mit Umgebungsvariablen

os.linesep # Zeilenseparator; \r\n win , \n linux

os.name # Name des OS -spez. Moduls (posix , dos ,...)

• Some general purpose methods

os.getcwd () # current working directory

os.chdir(path) # wechseln ins Verzeichnis path

os.getgroups () # Gruppen (UNIX)

os.uname() # Systeminformation (UNIX)

os.times() # (usr , sys , c_usr , c_sys , wct)

...




OS Methoden fur Dateien

chmod(path , mode) # chmod

chmod("tmp.py", 0664) # Modus erfordert vier Zeichen

listdir(path) # ls

mkdir(path) # mkdir

makedirs(path) # mkdir -p

rename(src , dst) # mv

remove(path) # rm

rmdir(path) # rmdir

removedirs(path) # rm -rf

stat(path) # stats (atime , ...)




Beispiel: AktienkurseAufgabe:• In mehreren csv-Dateien liegen Kursdaten vor• Eine Datei pro Firma (Dateiname=Firmenname)• Spalte mit den Kursen wird ermittelt• Methode zum Einlesen der Daten in interne Datenstruktur• Grafische Darstellung der Daten• Ermittlung von Mittelwerten

Folgenes wird geubt:• Kontrollstrukturen• Funktionen• Module• Dictionaries• Dateizugriff• Stringmanipulation mit s.split(chr)• Reduce-Operationen (min, max, sum)




Format der csv-Dateien (finance.yahoo.com)

Date ,Open ,High ,Low ,Close ,Volume ,Adj Close

2010 -11 -01 ,66.68 ,71.43 ,65.55 ,68.10 ,90700 ,68.10

2010 -10 -01 ,63.15 ,70.10 ,59.55 ,66.00 ,304000 ,66.00

2010 -09 -01 ,50.53 ,64.65 ,50.38 ,63.30 ,120700 ,63.30

2010 -08 -02 ,53.80 ,55.85 ,47.65 ,48.20 ,177600 ,48.20

2010 -07 -01 ,51.45 ,56.89 ,50.05 ,54.05 ,186700 ,54.05

2010 -06 -01 ,49.01 ,55.10 ,47.21 ,50.55 ,415500 ,50.55

2010 -05 -03 ,50.90 ,52.27 ,42.22 ,49.19 ,1485400 ,49.19

2010 -04 -01 ,47.58 ,53.01 ,46.60 ,50.94 ,1098400 ,50.94

2010 -03 -01 ,41.95 ,47.57 ,41.80 ,47.01 ,788800 ,47.01

2010 -02 -01 ,46.47 ,48.74 ,40.58 ,41.81 ,647900 ,41.81

2010 -01 -04 ,53.62 ,54.38 ,45.00 ,45.81 ,553800 ,45.81

2009 -12 -01 ,52.51 ,54.63 ,50.74 ,53.30 ,326000 ,53.30

2009 -11 -02 ,48.27 ,54.10 ,46.03 ,50.93 ,406200 ,50.93

2009 -10 -01 ,50.01 ,56.57 ,47.00 ,48.23 ,536200 ,48.23

...




Notige Module• os fur listdir, pylab zum plotten

import os

from pylab import *

Ermitteln aller Firmen• Im Verzeichnis liegen beliebig viele csv-Datein

def liesFirmen(dir):

firmenliste = []

dateien = os.listdir(dir)

for dateiname in dateien:

if dateiname [-4:] == ".csv":

firmenliste.append(dateiname [: -4])

return firmenliste




Liste mit Kurswerten einer Firma einlesen• Datei zum Lesen offnen• Datei zeilenweise durchgehen• Elemente der letzten Spalte extrahieren

def liesKurse(firma):

datei = open(firma+’.csv’, ’r’)

datei.readline () # Kopfzeile

preise = []

for zeile in datei:

spalten = zeile.split(’,’)

preis = spalten [-1]

preise.append(float(preis ))

datei.close()

preise.reverse ()

return preise




Automatisch alles einlesen• Firmennamen herausfinden• Dateien Einlesen• Speichern der Kurswerte in einem Dictionary• Berechnen von Mittelwerten

def init ():

daten =

firmen = liesFirmen(".")

for f in firmen:

daten[f] =

daten[f][’preise ’] = liesKurse(f)

daten[f][’avgpreis ’] = sum(daten[f][’preise ’]) \

/ len(daten[f][’preise ’])

daten[f][’minpreis ’] = min(daten[f][’preise ’])

daten[f][’maxpreis ’] = max(daten[f][’preise ’])

return daten




Plotten der Daten• Details zum Plotten gibt es spater

def zeichneKurse(daten , firmenliste ):

for firma in firmenliste:

plot(daten[firma ][’preise ’])

show()

• Verwendung als Programm und Modul

if __name__ == ’__main__ ’:

daten = init()

maxlen = max([len(key) for key in daten.keys ()])

for firma in daten.keys ():

name = firma+" "*(15-len(firma ))

min = daten[firma][’minpreis ’]

max = daten[firma][’maxpreis ’]

avg = daten[firma][’avgpreis ’]

print "%s: Preis(min/max/avg): %6g, %6g, %6g" \

%(name , min , max , avg)




• Ausfuhrung als Programm

>python aktienkurse.py

Deutsche_Bank (min/max/avg): 25.56 , 153.55 , 80.5186

Microsoft (min/max/avg): 15.66 , 42.8, 24.5753

Daimler (min/max/avg): 22.62, 110.15 , 47.6528

• Ausfuhrung als Modul

>>> from aktienkurse import *

>>> daten = init()

>>> zeichneKurse(daten , daten.keys ())




Teil VI

ObjektorientierteProgrammierung




Was ist ein Objekt?

• Ein Auto• Eine Katze• Ein Stuhl• ...

• Ein Molekul• Ein Stern• Eine Galaxie• ...

• Alles, was ein reales Objekt reprasentiert• Alles, was ein abstraktes Objekt reprasentiert• Alles, was irgendwelche Eigenschaften besitzt




Wie programmiert man objektorientiert?

• Wir mussen unsere Denkweise andern!• Wir mussen unseren Blickpunkt andern!• Es ist keine rein technische Frage.

Besonderheiten in python?• Wir verwenden schon die ganze Zeit objekorientierte Konzepte.• Wir merken es nur nicht.• In Python ist alles ein Objekt (sogar int).• Woher weiß z.B. print, wie ein int als Text aussieht?

>>> a=4

>>> a.__str__ ()

’4’




Theorie: Klassen und Objekte

Klassen:”Idee“ aller Objekte einer Art;Beispiel: Auto

Instanzen / Objekte:Konkretes Ding;Beispiel: ”Das rote Auto dort“




Theorie: Klassen und ObjekteKlassen: ”Idee“ aller Objekte einer Art; Bsp.: Auto• Eine Klasse ist so etwas ahnliches wie die Definition eines Typs.• Enthalt spezifische Eigenschaften der zu erzeugenden Objekte• Daten, die einmalig fur die Gesamtheit an Objekten gespeichert

werden (Klassenvariable / statische Variable)• Daten, die fur jedes (in jedem) Objekt gespeichert werden

(Objektvariable)• Operationen, die auf Objekten ausgefuhrt werden (Methoden)• Reprasentiert eine Klasse von Dingen/Objekten

Instanzen / Objekte: Konkretes Ding; Bsp.: ”Das rote Auto dort“• Eine Instanz ist ein konkretes Objekt, das zu einer Klasse gehort.• Eine Klasse spezifiziert Eigenschaften, ein Objekt hat konkrete

Werte fur diese Eigenschaften.• Zu einer Klasse konnen beliebig viele Objekte gehoren.• Jedes Objekt gehort genau zu einer Klasse (Ausnahme:

Polymorphismus)T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Grundstruktur einer Klasse

class ClassName:

<statement -1>

.

.

.

<statement -N>

Typische Statements• Funktionsdefinitionen⇒ methods• Variablendefinitionen⇒ data attributes (data members, fields, ...)




Theorie: Konstruktor, DestruktorBeim Erzeugen einer Instanz wird...• Speicher allokiert• eine Variable zum Referenzieren des Objekts erzeugt• vielleicht manches initialisiert• die Initialisierung mit einem Konstruktor (__init__(self,...))

durchgefuhrtWenn ein Objekt nicht mehr gebraucht wird, ...• wird Speicher freigegeben• mussen einige Dinge aufgeraumt werden• wird der Destruktor (__del__(self,...)) aufgerufen• ABER: Es ist nicht garantiert, wann er aufgerufen wird

class Example:

def __init__(self , ...):

do_something

def __del__(self , ...):

do_something




Theorie: self

Verwenden eines Objekts• Zugriff auf Methoden: objekt.methode(...)• Beispiel: liste.append(’x’)• Zugriff auf Variablen: objekt.variable• Beispiel: complex.real

Zugriff innerhalb der Objektmethoden• Konkretes Objekt wird außerhalb der Klassendefinition erzeugt⇒ Den Klassenmethoden ist der Name nicht bekannt⇒ Alle Methoden bekommen zusatzlichen Formalparameter self• self ist immer der erste Formalparameter, dieser wird beim

Aufruf (Aktualparameter) weggelassen• Innerhalb der Methoden Zugriff auf andere Methoden und

Variablen mit self.variable bzw. self.methode()




Erstes Klassenbeispiel

class Student:

"Einfache Beschreibung eines Studierenden"

def __init__(self , name): # Konstruktor

self.name = name

def getName(self): # Methode

return self.name

person1 = Student("Alex")

person2 = Student("Michaela")

print "Hier kommt " + person1.getName ()




Theorie: Methoden

• . . . implementieren das Verhalten von Objekten.• . . . reprasentieren alle Arten von Veranderungen eines Objekts

nach der Initialisierung.• accessor methods: geben Info uber Zustand des Objekts• mutator methods: verandern Zustand des Objekts⇒ verandern (mind.) eine Objektvariable




Theorie: Objektubergreifende (statische) Daten

Statische Variablen (Klassenvariablen)• Variablen pro Objekt speicherbar: Objektvariable

Beispiel: Name in Student

⇒ in Konstruktor definieren mit self• Variablen pro Klasse speicherbar: Klassenvariable⇒ nur 1x pro Klasse existentKlassisches Beispiel: Counter

Statische Methoden• Methoden, die lediglich statische Variablen (Klassenvariablen)

benutzen und unabhangig von konkreten Objekten verwendbarsein sollen

• Definition via Schlusselwort @staticmethod




Beispiel: Student reloaded

class Student:

"Einfache Beschreibung eines Studierenden"

anzahl = 0 # Klassenvariable

def __init__(self , name): # Konstruktor

self.name = name # Objektvariable

Student.anzahl += 1

def getName(self): # Methode

return self.name

@staticmethod

def getAnzahl (): # statische Methode

return Student.anzahl

person1 = Student("Alex")

person2 = Student("Michaela")

print "Hier kommt " + person1.getName ()

print "Anzahl Studis=", Student.getAnzahl ()




Prozeduraler Ansatz: Ofen

#!/usr/bin/python

def backe(temperature , mode):

Anweisungsblock

ofentemperatur = 180

ofenmodus = 2

backe(ofentemperatur , ofenmodus)

Bei einer ganzen Backerei?

ofentemperaturen = []

ofenmodi = []

ofentemperaturen.append (180)

ofenmodi.append (2)

...

for i in range(len(ofentemperaturen )):

backe(ofentemperaturen[i], ofenmodi[i])




Objektorientierter Ansatz: Ofen

#!/usr/bin/python

class Ofen(object ):

def __init__(self , temperatur , modus):

self.temperatur=temperatur

self.modus=modus

def backe(self):

Anweisungsblock

meinOfen = Ofen (180 ,2)

meinOfen.backe ()

Und wieder die ganze Backerei

ofenliste = []

ofenliste.append(Ofen (180 ,2))

...

for ofen in ofenliste:

ofen.backe ()




Was sind die Unterschiede?Prozeduraler Ansatz• Prozeduren/Funktionen legen fest, wie etwas ”produziert“wird• Bei der Implementierung denkt man an durchzufuhrende

Aktionen• Um die Speicherung der Daten muss man sich selbst kummern

Objektorientierter Ansatz• Objekte spezifizieren ”Dinge“(real oder abstrakt)• Bei der Implementierung muss man sich das Ding und seine

Eigenschaften vorhalten• Samtliche Daten sind im Objekt selbst gespeichert (und damit

gekapselt!)⇒ Wissen uber das Was, aber nicht uber das Wie⇒ Modularisierung! (großere) Unabhangigkeit versch. Code-Teile• Wenn viele Daten/Variablen mit einer Aktion verbunden sind,

lohnt sich der objektorientierte Ansatz besonders




Informationskapselung / Information Hiding

public• per default erstmal alles (alle Variablen/Methoden)• komfortabel (Zugriff von außen)• oft gefahrlich / unerwunscht

private• Zugriff von außen (d.h. außerhalb der Klasse) verboten⇒ Variablen/Methoden außen nicht sichtbar⇒ Variablen/Methoden außen nicht nutzbar• Workaround in python: name mangling⇒ Syntax: beginnend mit __, aber nicht endend mit __

(__privateVar, __privateMet())• Achtung: indirekt zugreifbar (._Klasse__privateVar)




Spezielle Variablen und Methoden

• __get__(value): Erlaubt Lesezugriff mit []• __set__(i, value): Erlaubt Schreibzugriff mit []• __add__(value): Erlaubt addition mit +• __str__(): Wird von print verwendet• __call__(): Erlaubt Aufruf mit objekt()

class Test:

def __init__(self , val):

self.value = val

def __add__(self , b):

return self.value + b.value

def __mul__(self ,b):

return self.value * b.value




Vererbung in der Realitat

Objekt Objekt

Maschine

Tier

Möbel

Pflanze

Comput. Objekt

Maschine

Auto

Tier

Säuget.

Möbel

StuhlTisch

Pflanze

Baum

Blume

Comput.

Laptop

Desktop

Vectorc.

Cluster

Reptilien

Objekt

Maschine

Auto

BMW

Audi

VW

Tier

Säuget.

Pferd

Esel

Möbel

StuhlTisch

Pflanze

Baum

Blume

Comput.

Laptop

Desktop

Vectorc.

Cluster

Liegestuhl

Ahorn

Eiche

Buche

Lilie

Tulpe

Rose

Reptilien

Katze

Krokodil

Schlange

Objekt

Maschine

Auto

BMW

Audi

VW

Golf

Polo

Phaeton

Tier

Säuget.

Pferd

EselMaulesel

Möbel

StuhlTisch

Pflanze

Baum

Blume

Comput.

Laptop

Desktop

Vectorc.

Cluster

Liegestuhl

Ahorn

Eiche

Buche

Lilie

Tulpe

Rose

Reptilien

Katze

Krokodil

Schlange




Vererbung• Klassen konnen von einer Basisklasse / Oberklasse erben• Dadurch sind sie abgeleitete Klasse / Unterklasse der

Basisklasse• Alle Klassen konnen spezialisiert werden• Alle Klassen zusammen bilden eine Klassenhierarchie• Methoden der Oberklasse konnen ubernommen oder neu

definiert (uberschrieben) werden

class Werkstudent(Student ):

def __init__(self , name , thema):

Student.__init__(self , name)

self.thema = thema

def printInfo(self):

print "%s arbeitet an: %s" % (self.name , self.thema)

person3 = Werkstudent("Stefan", "Tabellenkalkulation")

person3.printInfo ()




Modellierung mit UML

Flussdiagramme• Nur geeignet fur kleine Algorithmen• Keine Darstellung von Klassenhierarchien• Keine Darstellung von Datenabhangigkeiten• ...

UML: Unified Modeling Language•

”visuelle“Modellierungssprache• Sprache unterstutzt unterschiedliche Diagrammtypen• Software wird uber diese Diagramme/Modelle spezifiziert• Kann als Basis fur die Dokumentation dienen• Automatische Code-Generierung moglich• Auch fur sonstige betriebliche Ablaufe geeignet




Strukturdiagramme• Klassendiagramm• Objektdiagramm• Kompositionsstrukturdiagramm• Komponentendiagramm• Verteilungsdiagramm• Paketdiagramm• Profildiagramm

Verhaltensdiagramme• Aktivitatsdiagramm• Anwendungsfalldiagramm (Use-Cases)• Interaktionsubersichtsdiagramm• Kommunikationsdiagramm• Sequenzdiagramm• Zeitverlaufsdiagramm• Zustandsdiagramm




KlassendiagrammEinzelne Klasse• Klasse dargestellt durch Rechteck• Horizontale Unterteilung in drei Bereiche• Oberster Bereich: Klassenname• Mittlerer Bereich: Attribute/Variablen• Unterer Bereich: Methoden

Klassenname

+ public_vars private_vars+ public_methods private_methods




Syntax der Klassenbeschreibung• Name von Klasse/Variablen/Methoden frei wahlbar• +: public Variable/Methode (optional)• -: private Variable/Methode (optional)• Fur Variablen kann der Typ angegeben werden (nach :)• Fur Methoden kann der Ruckgabewert angegeben werden (nach

:)




Vererbung• Darstellung durch Pfeil mit ”geschlossener“Spitze•

”ist ein“- Beziehung• Pfeil von spezialisierter Klasse zur Vaterklasse

Fahrzeug

Auto Motorrad




Assoziation• Semantischer Zusammenhang zwischen Klassen• Beschreibung der Assoziation neben dem Pfeil

Vorlesung

Dozent Student

liest hört




Aggregation und Komposition• Aggregation (leere Raute): ”hat ein“• Komposition (volle Raute): ”enthalt ein“• Angabe der Anzahl moglich

Fahrzeug

Besitzer Türen1 5




Teil VII

Regulare Ausdrucke




Was wir uns (vielleicht) schon immer gefragt haben:• Wie funktioniert Suche nach einer Zeichenfolge in einem Text?• Wie wird Auto-Vervollstandigung in Entwicklungsumgebungen

realisiert?• Woher weiß ein email client, ob eine Emailadresse valide ist?

⇒ Da sind regulare Ausdrucke im Spiel

regulare Ausdrucke = regular expressions = regex

Unterbau: Formale Sprachen (Theoretische Informatik)




Formale SprachenSyntax• Jede Sprache (naturliche und formale) hat Regeln.• Fur einen Text lasst sich feststellen, ob er zu einer Sprache

gehort:• Deutsch: Die Studierenden sitzen in der Vorlesung.• Python: print "Hello World"

• Eine Grammatik beschreibt, nach welchen Regeln eine Sprachesyntaktisch aufgebaut ist.

• Bei einer Programmiersprache pruft der Compiler/Interpreter dieSyntax.

Semantik• Auch bei grammatikalisch korrektem Aufbau kann ein Text

sinnlos sein• Deutsch: Die Vorlesung sitzt in den Studierenden.• Python: print "sin(x) = %f"% cos(x)

• Eine Semantik-Korrektur fur Programmiersprachen gibt es nochnicht.




Formale Grammatik

• Eine formale Grammatik ist ein 4-Tupel G = (V ,Σ,P,S)

• V ist eine endliche Menge, das Vokabular• Σ ⊂ V ist das Alphabet, die Elemente heißen Terminalsymbole

(”Student“, ”Vorlesung“; a, b, ...)• N = V\Σ sind die Nichtterminalsymbole/Variablen (<Subjekt>,<Verb>; A, B, ...)

• X ∗ ist die Kleensche Hulle der Menge X . Enthalt beliebigeKonkatenationen von Elementen aus der Menge.

• P ⊂ (V ∗\Σ∗)xV ∗ ist die Menge der Produktionsregeln.Uberfuhrung eines Wortes/Texts R, das mindestens einNichtterminal enthalt (R ∈ V ∗\Σ∗) in ein beliebiges Wort Q ∈ V ∗

<Subj ><Verb ><Obj > --> Der Student <Verb ><Obj >

AB --> AaB

• S ∈ (V\Σ) ist das Startsymbol




Chomsky-Hierarchie

• Eingefuhrt von Noam Chomsky• Eine Hierarchie von Klassen formaler Grammatiken

Typ 0: unbeschrankte Grammatik• enthalt alle formalen Grammatiken• zugehorige Sprache wird von Turingmaschine akzeptiert

Typ 1: kontextsensitive Grammatik• Produktionsregeln der Form αBγ → αβγ (B Nichtterminal,

griechische Buchstaben Worte aus V ∗)• Sprache wird von linear beschrankter Turingmaschine erkannt

Typ 2: kontextfreie Grammatik• Produktionsregeln der Form A→ α

• erkannt von Kellerautomat, Programmiersprachen sind Typ 2Typ 3: regulare Grammatik• Produktionsregeln der Form A→ a und A→ aB• erkannt durch endliche Automaten / regulare Ausdrucke




Regulare Grammatik der Exponentialzahlen• Einschrankende Annahmen:

• Beide Vorzeichen (Anfang+Exponent) zwingend anzugeben• Genau eine Vor- und mindestens eine Nachkommastelle

• G = (V ,Σ,P,S)• V = Σ ∪ N• N = S,M,N1,P,N2,E1,E2• Σ = −,+,0,1,2,3,4,5,6,7,8,9, .,e• Platzhalter: z ∈ 0− 9

Produktionsregeln (rechts Beispiel: -3.81e-11)

S -> -M S -> +M | S -> -M

M -> zP | -> -3P

P -> .N1 | -> -3.N1

N1 -> zN2 | -> -3.8N2

N2 -> zN2 N2 -> eE1 | -> -3.81N2 -> -3.81eE1

E1 -> -E2 E1 -> +E2 | -> -3.81e-E2

E2 -> zE2 E2 -> z | -> -3.81e-1E2 -> -3.81e-11




Regulare Sprachen und Ausdrucke???Die Menge der regularen Sprachen uber einem Alphabet Σ und diezugehorigen regularen Ausdrucke sind rekursiv definiert:• Die leere Sprache Ø ist eine regulare Sprache und der

zugehorige regulare Ausdruck ist Ø.• Der leere String ∧ ist eine regulare Sprache und der

zugehorige regulare Ausdruck ist ∧.• Fur jedes a in Σ, ist die einelementige Sprache a eine

regulare Sprache und der zugehorige regulare Ausdruck ist a• Wenn A und B regulare Sprachen sind mit zugehorigen

regularen Audrucken r1 and r2, dann• ist A ∪ B (Vereinigung) eine regulare Sprache und der

zugehorige regulare Ausdruck ist (r1|r2)• ist AB (Verknupfung) eine regulare Sprache und der

zugehorige regulare Ausdruck ist (r1r2)• ist A∗ (Kleene star) eine regulare Sprache und der

zugehorige regulare Ausdruck ist (r1∗)




Regularer Ausdruck fur Exponentialzahlen• Es gibt unendlich viele Sprachen uber jedem endlichen

nichtleeren Alphabet• z.B. fur Σ = a sind Sprachen: , a, aa, a,aa, aaa, ...• Wir wollen zeigen, dass sich die Sprache der Exponentialzahlen

und der zugehorige regulare Ausdruck rekursiv herleiten lassen• Alphabet: −,+,0,1,2,3,4,5,6,7,8,9, .,e• Einelementige Sprachen: S− = −, S+ = +, S0 = 0,

S1 = 1, ..., S. = ., Se = e• zugehorige regularen Ausdrucke:

r− = −, r+ = +, r0 = 0, r1 = 1, ..., r. = ., re = e• Vereinigung zweier Sprachen:a,b ∪ 0,1,2 = a,b,0,1,2 ⇒ |S1|+ |S2| Elemente

• Verknupfung zweier Sprachen:a,b0,1,2 = a0,a1,a2,b0,b1,b2 ⇒ |S1| · |S2| Elemente

• Kleensche Hulle einer Sprache:0,1∗ = Ø,0,1,00,01,10,11,000, ... ⇒ ∞ Elemente




Erzeugung weiterer Sprachen aus den einelementigenOperation Sprache AusdruckS− ∪ S+ SA = −,+ rA = (−|+)S0 ∪ S1∪ SB = 0, 1, ..., 9 rB = (0|1|...|9) = (0− 9)SASBS. SC = −0., ...,−9., ... rC = (−|+)(0− 9).SCSB SD = −0.0,−0.1, ... rD = (−|+)(0− 9).(0− 9)SDS∗

B SE = −0.0,−0.00, ... rE = (−|+)(0− 9).(0− 9)(0− 9)∗

rE = (−|+)(0− 9).(0− 9)+

SE Se SF = −0.0e, ... rF = (−|+)(0− 9).(0− 9)+eSF SA SG = −0.0e−, ... rG = (−|+)(0− 9).(0− 9)+e(−|+)SF SB SH = −0.0e − 0, ... rH = (−|+)(0− 9).(0− 9)+e(−|+)(0− 9)SHS∗

B SI = −0.0e − 00, ... rI = (−|+)(0− 9).(0− 9)+e(−|+)(0− 9)+

• SI ist genau die Sprache der Exponentialzahlen, fur die bereitseine Grammatik erstellt wurde

• rI = (−|+)(0− 9).(0− 9)+e(−|+)(0− 9)+ ist der zugehorigeregulare Ausdruck

• in Python: rI = r’[-+][0-9]\.[0-9]+e[-+][0-9]+’

• noch kurzer: rI = r’[-+]\d\.\d+e[-+]\d+’




Regulare Ausdrucke!!!What does all this mean to you, as a user? Absolutely nothing. As auser, you don’t care if it’s regular, nonregular, unregular, irregular, orincontinent. So long as you know what you can expect from it, youknow all you need to care about. — Jeffrey FriedlWo fast jeder sie schon verwendet hat:• Suche nach einer Zeichenfolge in einem Text• Tabulator-Vervollstandigung• Mit * eine Gruppe von Dateien auszuwahlen (*.txt)

Einschrankungen• Fur Anfanger sehr kryptisch• Nach reiner Lehre ”Zahlen“ nicht moglich (z.B. anbn geht nicht)

Vorteile• Sehr nutzliches Werkzeug zum Finden von Pattern• Vergleichbare ”manuelle“ Implementierung viel aufwandiger• Python hat Erweiterungen, die sogar ”Zahlen“ ermoglichen




Regulare Ausdrucke in PythonPattern Syntax• Regulare Ausdrucke immer als ”raw string“• E.g. r’([0-9]*\.[0-9]*|[0-9]*)’. entspricht beliebigem Zeichen außer newline^ entspricht dem Beginn eines strings$ entspricht dem Ende eines strings* voriger Ausdruck beliebig oft (incl. Null mal) (gierig)+ voriger Ausdruck beliebig oft (excl. Null mal) (gierig)? voriger Ausdruck Null- oder einmal (gierig)*?, +?, ?? nicht gierige Versionen von *, +, ?

m voriger Ausdruck genau m Malm,n voriger Ausdruck m bis n Mal (gierig)m,n? nicht gierige Version von m,n[...] ein beliebiges Zeichen aus der Menge (...)[^...] ein beliebiges Zeichen, das nicht in der Menge istA|B A oder B (A und B sind regulare Ausdrucke)(...) speichert den gefundenen Inhalt der Klammern




Regulare Ausdrucke - Maskierungszeichen

• Zeichen mit einer speziellen Bedeutung (., *, ...) konnen mitihrer eigentlichen Bedeutung durch voranstellen eines\verwendet werden, z.B. \., \*

• Normale Maskierungszeichen (\n, \t, ...) funktionieren wieerwartet, z.B. r’\n+’ entspricht einem oder mehrerenZeilenumbruchen

\number enspricht dem n-ten zuvor gefundenen Text (starting from 1)\d entspricht [0-9]\D entspricht [^0-9]\s enspricht beliebigem whitespace (= [\t\n\r\f\v])\S alles außer whitespace\w beliebiges alphanumerisches Zeichen (= [a-zA-Z0-9 ])\W beliebiges nicht-alphanumerisches Zeichen\A entspricht dem Beginn eines strings\Z entspricht dem Ende eines strings




Regulare Ausdrucke - erste Python-Beispiele

• findall(patt, string) – finde alle nicht uberlappendenVorkommen von patt in string

• sub(patt, repl, string) – ersetze Vorkommen durch repl

import re

regstr = r’\d+\.\d+|\d+’

s = "12.3/5/4.4;5.7;6"

re.findall(regstr ,s)

regstr = r’\d*\.\d*|\d*’

s = "12.3/5/4.4;5.7;6"

re.findall(regstr ,s)

regstr = r’(\d*\.\d*|\d*)’

re.sub(regstr , r’\1xxx’, s)




Regulare Ausdrucke - gierig/greedy

• Was bedeutet gierig/nicht-gierig (bzw. greedy/non-greedy)?• Beispiel:

reg_greedy = r’<.*>’ #match as many chars as possible

s = "<H1 >title </H1 >"

findall(reg_greedy ,s) #whole ’<H1 >title </H1 >’ matched

reg_nongreedy = r’ <.*?>’ #match as few chars as possible

findall(reg_nongreedy ,s) #only ’<H1 >’ matched

• gierig/greedy: so viele Zeichen wie moglich gematched; nur,wenn das schiefgehen wurde, gibt es ein backtracking durch dieregex engine




Regulare Ausdrucke - Match-Objekte

• Manche Funktionen aus re geben keine strings, sondern

”Match-Objekte“(m) zuruck• match(patt, string) Sucht nach Ubereinstimmung am Beginn

des strings• search(patt, string) Sucht die erste Ubereinstimmung im string• finditer(patt, string) wie findall, gibt aber einen iterator zuruck• m.group(i) Gibt Ubereinstimmungsgruppe i zuruck (i ≥ 1)• m.groups() Gibt Tupel mit allen Ubereinstimmungsgruppen

zuruck

import re

s = "Ferien 12/24/2010 - 01/06/2011"

patt = r’(\d+)/(\d+)/(\d+)’

for m in re.finditer(patt ,s):

print("%s.%s.%s" % (m.group(2), m.group(1), m.group (3)))




Regulare Ausdrucke - Ratsel zu Primzahlen

Wieso “berechnet” folgender Code Primzahlen?

import re

def isPrime(p):

m = re.search(r’^1?$|^(11+?)\1+$’, p)

if(m):

return False

else:

return True

for i in xrange (0 ,1024):

if(isPrime(’1’ * i)):

print i




Teil VIII

Exceptions




Exceptions: Fehler abfangen

• Fehler treten beispielsweise auf, wenn ein Benutzer interagiert:

x = raw_input("Bitte eine Zahl eingeben: ")

print float(x)

• Wenn die Konvertierung nach float nicht moglich ist, wird eineValueError exception geworfen

Bitte eine Zahl eingeben: xyz

...

Traceback (most recent call last):

File "<stdin >", line 1, in <module >

ValueError: invalid literal for float (): xyz

• Um mit diesem fehlerhaften Verhalten umzugehen, kann dieexception abgefangen und bearbeitet werden




Syntax fur Exceptions• try: leitet einen Block ein, der normal ausgefuhrt wird.• Darauf folgen (mehrere) except-Anweisungen, die im

ensprechenden Fehlerfall ausgefuhrt werden.• Dann optional ein else, dass nur ausgefuhrt wird, wenn kein

Fehler auftritt;• Zuletzt optional finally, das immer ausgefuhrt wird.

try:

# beliebiger Code

except ValueError as e: # Py >3.0: "Exception as e"

# behandle ValueError

except Exception as e:

# behandle alle ubrigen Fehler

else:

# ausfuhren , falls keine Fehler auftrat

finally:

# Wird immer ausgefuhrt




BeispielWieder das Beispiel mit User-Eingabe

eingabe = False

while not eingabe:

x = raw_input("Bitte eine Zahl eingeben: ")

try:

print float(x)

eingabe = True

except ValueError as e:

print "Das war keine Zahl!"

# evtl. noch irgendwas mit e machen

• Eigene Exceptions werfen (Hierfur konnen Standardfehlerverwendet werden, oder eigene Fehlertypen erzeugt werden):

raise RuntimeError("Ooops! Irgendwas ging schief!")

raise IOError("Datei existiert nicht ...")




Eingebaute Exceptions

BaseException # Wurzel aller exc.

GeneratorExit

KeyboardInterrupt # z.B. Ctrl+C

SystemExit # Program Exit (sys.exit ())

Exception # Basis fur nicht -exit exc.

StopIteration

StandardError # Nur Python 2.x

ArithmeticError

FloatingPointError

ZeroDevisionError

AssertionError

AttributeError

EnvironmentError

IOError # z.B. bei open(" datei.txt")

OSError

EOFError

ImportError # z.B. Modul nicht gefunden




LookupError

IndexError # z.B. bei Tupeln/Listen

KeyError # z.B. bei Dictionaries

MemoryError

NameError # Name nicht gefunden

UnboundedLocalError

ReferenceError

RuntimeError

NotImplementedError

SyntaxError

IndentationError

TabError

SystemError

TypeError # Typ -Fehler (2 + "3")

ValueError # Wert -Fehler (float("drei "))

UnicodeError

UnicodeDecodeError

UnicodeEncodeError

UnicodeTranslateError




Built-In Exception werfen

def readfloat1 ():

while True:

try:

a = raw_input("Zahl zwischen 0.0 und 1.0: ")

a = float(a)

if(a<0.0 or a >1.0):

raise ValueError("Zahl nicht in [0.0; 1.0]")


print "Fehler: %s" % e

else:

return a

• Bei a = float(a) wird automatisch Fehler geworfen, der Rest destry-Blocks wird nicht mehr ausgefuhrt

• Falls a konvertierbar ist, wird der Wertebereich uberpruft und ggfneue Exception geworfen




Verschachtelte Try-Blocke

def readfloat2 ():

while True:

try:

a = raw_input("Zahl zwischen 0.0 und 1.0: ")

try:

a = float(a)


raise ValueError("Das ist keine Zahl: %s" % a)

if(a<0.0 or a >1.0):

raise ValueError("Zahl nicht in [0.0; 1.0]")


print "Fehler: %s" % e

else:

return a

• Exceptions konnen abgefangen und gleich wieder geworfenwerden.




Eigene Exceptions

• Hierfur sind Klassen notig (vgl. Kap. 6)• Kleiner Ausblick:

class MeinFehler(Exception ): pass

raise MeinFehler("schlimmer Fehler")

• Man kann naturlich auch kompiziertere Dinge tun:

class KomplizierterFehler(Exception ):

def __init__(self , nummer , nachricht ):

self.args = (nummer , nachricht)

self.nummer = nummer

self.nachricht = nachricht

def machIrgendwas ():

...




Teil IX

2D-Grafiken mit TKInter




Visualisierung

R.W. Hamming:”The purpose of computing is insight, not numbers.“




Rekursion• Bei einer Rekursion wird eine Funktion durch sich selbst

definiert.• Beim Aufruf der Funktion wird die Funktion selbst wieder

aufgerufen.• Abbruchkriterium verhindert endlose Rekursion

Fibonacci-Folge (1202)• Fibonacci beschrieb damit das Wachstum von

Kaninchenpopulationen:• fn = fn−1 + fn−2• f0 = 0• f1 = 1

def fib(n):

if n<=1:

return n

else:

return fib(n-1)+ fib(n-2)




Koch-Kurve (1904)

• Eines der ersten entdeckten Fraktalen Objekte• Uberall stetig, nirgends differenzierbar• Beginn mit einer Strecke• Mittleres Drittel wird entfernt• Stattdessen zwei gleichlange Strecken, die ein Dreieck bilden

• Kann ebenfalls mit Rekursion erzeugt werden• Wir benotigen eine Funktion, die die Kurve malt.• Um die Kurve zu malen, mussen wir erst alle vier Teilstrecken

malen.• Fur jede Teilstrecke, mussen deren vier Teilstrecken gemalt

werden.• ...




Visualisierung der Koch-Funktion mit turtle

import turtle

def kochkurve(laenge , ebene):

if(ebene >0):

kochkurve(laenge /3.0, ebene -1)

turtle.left (60)


turtle.right (120)


turtle.left (60)


else:

turtle.forward(laenge)




Initialisierung

def init ():

turtle.reset ()

turtle.setworldcoordinates (-0.1,-0.9, 1.1, 0.3)

turtle.speed (0)

turtle.hideturtle ()

Ausfuhrung

if __name__ == ’__main__ ’:

n=int(raw_input("Anzahl Rekursionen:"))

init()

kochkurve (1.0, n)

turtle.right (120)

kochkurve (1.0, n)

turtle.right (120)

kochkurve (1.0, n)

turtle.exitonclick ()




Randnotiz: Raumfullende Kurven

• Analoge Erzeugung uber Rekursion• Zusatzlich: Grammatik notig (vgl. re)• Konkrete Nutzung: Parallelisierung, Cache-Nutzung• Beispiele: Hilbert-Kurve, Peano-Kurve




Koordinatenbasierte GUIs

GUIs fur Python• GUI: Graphical User Interface• Inzwischen verschiedene GUIs verfugbar

• TKinter: Python Schnittstelle fur Tcl/TK (im Standardenthalten)

• EasyGUI: baut auf TKinter auf, einfach aber beschrankt• PyGTK: basiert auf GTK+ (Linux)• wxPython: benutzt OS-Grafikbibliotheken• PythonCard: baut auf wxPython auf, einfach aber

beschrankt• PyQt: basiert auf Qt (z.B. KDE)• PythonWin: nur fur Windows• ...




Tcl/TK

• Tool Command Language / ToolKit• Relativ einfach• Ausreichend fur mittelgroße Anwendungen

widgets• Gangiges Wort fur jegliche GUI-Komponente• frame Kontainer fur alles mogliche, z.B. ein Fenster• canvas Leinwand/Zeichenflache• button

• checkbox

• radiobutton

• label z.B. Text• menu

• scrollbar

• ...




Module Tkinter

Hello world

import Tkinter as tk

root = tk.Tk()

root.title("Hello -Fenster")

label = tk.Label(root , text="Hallo Welt",

font="times 32 bold")

label.pack()

button = tk.Button(root)

button["text"] = "Hallo -Knopf"

button["background"] = "blue"

button.pack()

root.mainloop ()

• tk.Tk() erzeugt das Hauptfenster.• root.title(...) legt den Titel fest.• widget.pack() macht das widget sichtbar.• root.mainloop() wartet auf Ereignisse.




Prinzipielle Vorgehensweise• Festlegen, wie die GUI aussehen soll• Zugrundeliegende Funktionalitat implementieren• GUI-Elemente mit Funktionalitat verbinden

Hello world als Klasse


class MyVis:

def __init__(self , root):

self.label = tk.Label(root , text="Hallo Welt",

font="times 32 bold")

self.label.pack()

self.button1 = tk.Button(root)

self.button1["text"]= "Hallo -Knopf"

self.button1["background"] = "blue"

self.button1.pack()

root = tk.Tk()

myvis = MyVis(root)

root.mainloop ()




Anderung der Anordnung


class MyVis2:


self.button1 = tk.Button(root ,

text="Hallo -Knopf", background="blue")

self.button1.pack(side = tk.LEFT)

self.button2 = tk.Button(root)

self.button2.configure(text="Exit",

background="red")

self.button2.pack(side = tk.LEFT)

root = tk.Tk()

myvis = MyVis2(root)

root.mainloop ()

• Button-Optionen konnen bei der Konstruktion angegebenwerden.

• Button-Optionen konnen mit configure geandert werden.• side-Option fur pack bestimmt Ausrichtung




Aktionen aufrufen - Binding• Bestimmte Ereignisse erzeugen ein ”event“.• Z.B. Tastendruck, Mausklick, Mausbewegung, ...• Ereignisse konnen an widgets (Buttons,...) gebunden werden.• Der Knopf wartet dann auf das jeweilige Ereignis.• Tritt es auf, wird die assoziierte Funktion aufgerufen.

class MyVis2:


...

self.root = root

self.button1.bind("<Button -1>", self.action1)

self.button2.bind("<Button -1>", self.action2)

def action1(self , event):

print "Hallo Welt!"

self.button1["background"] = "red"

def action2(self , event):

self.root.destroy ()




Aktionen aufrufen - command• Bisher: Knopf reagiert nicht auf kompletten Klick• Event "<Button-1>" ist gleich "<ButtonPress-1>"

• Loslassen entspricht "<ButtonRelease-1>"• Mit command werden mehrere Events an ein widget gebunden

class MyVis3:



text="Hallo -Knopf", background="blue",

command=self.action1)


command=self.action2)

...

def action1(self):

...

def action2(self):

...




Parameter an Aktionen ubergeben• Bisher z.B.: command=self.action1• command ist die Funktion (deren Name)• Mit Parameter: command=self.action1(x)?• Hier ist command der Ruckgabewert der Funktion (None)• Losung: Funktion ohne Parameter angeben, die eine Funktion

mit Parametern aufruft (⇒ lambda-Funktion)

command=lambda: self.action1("Hallo Welt")

• Wert der Variablen zum Zeitpunkt des Events wird verwendet!

command=lambda: self.action1(var1 , var2)

• Variablenwerte zum Zeitpunkt der Erstellung:

command=lambda x=var1 , y=var2: self.action1(x, y)

• Funktionsheader: def action1(self, x, y):




widget Hierarchie

class MyApp:


# buttons frame mit einem Knopf

self.buttons_frame = tk.Frame(root)

self.buttons_frame.pack(side=tk.TOP)

self.button1 = tk.Button(self.buttons_frame)

self.button1.pack(side=tk.LEFT)

# top frame

self.top = tk.Frame(root)

self.top.pack(side=tk.TOP)

# linker + rechter Frame (Kinder von top)

self.left_frame = tk.Frame(self.top , background="blue",

borderwidth =5, height =200, width =80)

self.left_frame.pack(side=tk.LEFT) ###

self.right_frame = tk.Frame(self.top , background="red",

borderwidth =5, relief=tk.RIDGE , width = 220)

self.right_frame.pack(side=tk.RIGHT , fill=tk.BOTH)




widget-Hierarchie• Bisher wurden alle widgets direkt als Kinder von root erzeugt.• widgets konnen selbst Kinder haben.• Beliebige Verschachtelungen moglich• Fur alle widgets kann ein Padding (”Auffullen“) angegeben

werden.Geometrie-Manager: pack• Anordnung im Wesentlichen uber side=... (LEFT, TOP,...)

• fill gibt an, wie sich die Objekte ausdehnen sollen• tk.NONE gar nicht (minimale Ausdehnung)• tk.X, tk.Y nur in X/Y-Richtung• tk.BOTH in beide Richtungen

• expand Das widget versucht, moglichst viel Flache zu belegen• anchor Damit kann das widget sich in einem bestimmten Teil der

Flache platzieren (tk.N, tk.NW, tk.NE, ..., tk.CENTER)Geometrie-Manager: grid• Anordnung der widgets in einem Gitter• Z.B. widget.grid(row=2, column=7)




Zeichnen mit TkinterDas canvas Widget• Kann wie jedes widget in ein Frame gesetzt werden• Kann auch ohne Frame erzeugt werden:

cv = tk.Canvas(width =800, height =600)

cv.pack()

• Kann zum Zeichnen verwendet werden:

cv.create_line (50,50,650, 300)

cv.create_polygon (250, 50, 500, 500, 50, 500)

cv.create_text (400, 100, text="Hallo Welt")

• Weitere Methoden:• create_oval

• create_bitmap

• ...




Animation mit Tkinter

#! /usr/bin/python

from Tkinter import *

import time , math

cv = Canvas(width =800, height =600)

cv.pack()

oval = cv.create_oval (300, 50, 400, 150)

cv.itemconfig(oval , fill=’blue’)

i = 0.0

while True:

cv.move(oval , 10* math.cos(i), 10* math.sin(i))

cv.update ()

time.sleep (0.02)

i += 0.05




Teil X

Simulation von Partikelsystemen




Wo stehen wir?

Teil I: Erste SchritteTeil II: DatentypenTeil III: KontrollstrukturenTeil IV: Funktionen und ModuleTeil V: IO und Datentypen

Teil VI: Objektorientierte ProgrammierungTeil VII: Regulare AusdruckeTeil VIII: ExceptionsTeil IX: GrafikTeil X: PartikelsystemeTeil XI: Datenstrukturen

Teil XII: Wissenschaftliches RechnenTeil XIII: WarmeleitungTeil XIV: Losung von Linearen GleichungssystemenTeil XV: Mehrgitterverfahren




Simulationspipeline

Phänomen, (physikalischer) Vorgang

Mathematisches Modell (z.B. DGL)

Numerischer Algorithmus (z.B. Euler)

Simulations-Programm

Ergebnisse (Messwerte, Bilder, ...)

Modellbildung

Numerik (Diskretisierung, ...)

Implementierung

Visualisierung

Va

lidie

run

g




Simulationspipeline - reloadedTechnische Universität München

H.-J. Bungartz: The Bavarian Graduate School of Computational Engineering www.bgce.de BGCE – CSE – come.tum Joint Graduation Ceremony, November 10, 2011

2

The Simulation Pipeline Basis of Computational Engineering

Mathematical model Discretization & solver

Parallel implementation, HPC

Exploration

Validation

Software

Insight, Design




Simulation des Sonnensystems

• Phanomen: Planeten und Sterne ziehen sich an• Modell: Gravitationspotential/kraft fur jedes Paar:

Fi =∑

j 6=i Fij

• Newton:

ri =Fi

mi=

∑j 6=i Fij

mi= −

∑j 6=i

∂U(ri ,rj )∂|rij |

mi

• Numerik: Diskretisierung der DGL (ODE 2.Ordnung)⇒ Trafo auf System 1. Ordnung




Simulation des Sonnensystems II

Simulationsprogramm• Datenstrukturen fur die Objekte• Berechnung der Kraft• Losung der Bewegungsgleichung• Kopplung zur Visualisierung• Visualisierung mit Tkinter

Anmerkungen• Eine Vorlesung fur ein komplettes Programm ist wenig!• Das Programm wird funktional, aber weder effizient noch an

allen Stellen sehr sauber• In den Code-Stucken sind kaum Kommentare, diese werden

mundlich in der Vorlesung gegeben. Bei eigenen Programmengehoren die Kommentare in den Code!




UML-Diagram des Astro-Simulationsprogramms

Euler

+ integrate

Integrator- dt: float

+ integrate

Spherevis- dt: float

- pos2pixel+ drawSpheres

SimulationInter+ spheres+ canvas, dt

- calculateForces+ setIntegrator+ initVis+ simLoop GravSphere

+ calcForce

Sphere+ rad, mass, pos+ vel, force, color- dist+ resetForce+ calcForce

AstroInter

+ initVis

n

1

1




Integrator• Basisklasse Integrator als Schnittstellendefinition• Unterklasse Euler: Einfacher Integrator• Fur jedes Objekt Losung der Bewegungsgleichung in beiden

Raumrichtungen

class Integrator(object ):

def __init__(self , dt):

self.dt = dt

def integrate(self , sphereList ):

pass

class Euler(Integrator ):

def integrate(self , sphereList ):

for s in sphereList:

for d in range (2):

s.pos[d] = s.pos[d] + self.dt*s.vel[d]

s.vel[d] = s.vel[d] + self.dt*s.force[d]/s.mass




UML-Diagram (WH)

Euler

+ integrate


+ integrate





+ calcForce


AstroInter

+ initVis

n

1

1




Klasse zur Visualisierung: Spherevis

Konstruktor __init__

• Parameter fur den Konstruktor: canvas, posLL und posRR

• canvas enthalt die Zeichenflache• posLL und posUR enthalten Koordinaten der Ecken des

Simulationsgebiets• Konstruktor speichert ubergebenes canvas, dessen Große und

die Koordinaten

class Spherevis(object ):

def __init__(self , canvas , posLL , posUR):

self.canvas = canvas

self.width = canvas.winfo_reqwidth ()

self.height = canvas.winfo_reqheight ()

self.posLL = posLL

self.posUR = posUR




Methode pos2Pixel

• Zu zeichnende Objekte sollen in Koordinaten desSimulationsgebiets angegeben werden

• Diese mussen in Pixel umgerechnet werden• Dazu wird das Intervall des Simulationsgebiets auf das Intervall

zwischen Null und der Anzahl an Pixeln abgebildet• pos2Pixel erhalt als Parameter die Simulationskoordinaten und

gibt die Pixelkoordinaten zuruck

def pos2Pixel(self , posx , posy):

pix_x = self.width *(( posx - self.posLL [0])/

(self.posUR [0]-self.posLL [0]))

pix_y = self.height *(( posy - self.posLL [1])/

(self.posUR[1]-self.posLL [1]))

return pix_x , pix_y




Methode drawSpheres

• drawSpheres erhalt als Parameter eine Liste von Kugeln (KlasseSphere)

• drawSpheres loscht die Zeichenflache und malt alle Kugeln mitderen Radius und Farbe

• Fur jede Kugel wird die ”Bounding Box“in Pixeln ermittelt, dannmit create_oval der ensprechende Kreis gezeichnet

• Dann wird die Zeicheflache aktualisiert

def drawSpheres(self , spheres ):

self.canvas.delete("all")

for s in spheres:

pixelLL = self.pos2Pixel(s.pos[0]-s.rad ,

s.pos[1]-s.rad)

pixelUR = self.pos2Pixel(s.pos [0]+s.rad ,

s.pos [1]+s.rad)

self.canvas.create_oval(pixelLL [0], pixelLL [1],

pixelUR [0], pixelUR [1],

fill=s.color)

self.canvas.update ()




UML-Diagram (WH)

Euler

+ integrate


+ integrate





+ calcForce


AstroInter

+ initVis

n

1

1




Klassen fur Kugeln: Sphere und GravSphere


• Parameter: Position[3], Geschwindigkeit[3], Radius, Masse undFarbe

• Alle Parameter werden gespeichertMethode resetForce

• Setzt die Kraft (force) auf NullMethode dist

• Parameter: zweite Kugel• Gibt den Abstand (Vektor) zwischen den Kugeln zuruck

Methode calcForce

• Parameter: zweite Kugel• Berechnet die Kraft, die die zweite Kugel auf self ausubt und

addiert sie auf self.force• In der Basisklasse Sphere nur als ”Interface“• Konkrete Implementierung in Subklassen




Implementierung der Klasse Sphere

from math import sqrt

class Sphere(object ):

def __init__(self , pos , vel , rad , mass , color="red"):

self.rad = rad

self.mass = mass

self.pos = list(pos)

self.vel = list(vel)

self.force = [0,0]

self.color = color

def resetForce(self):

for d in range (2):

self.force[d] = 0

def dist(self , sphere2 ):

dist = [0,0]

for d in range (2):

dist[d] = sphere2.pos[d] - self.pos[d]

return dist

def calcForce(self , sphere2 ):

pass




Subklasse GravSphere• Gravitationskraft 1D: F (d) = G · m1m2

d2

• G: Gravitationskonstante (6.67428 · 10−11)• Gravitationskraft: ~F (~r) = G · m1m2

d3 ·~r• Dabei: ~r : Vektor zwischen beiden Kugeln, d : Abstand• GravSphere ubernimmt Konstruktor, resetForce und dist von der

Basisklasse• calcForce wird uberschrieben

class GravSphere(Sphere ):

gravConst = 6.67428e-11


dist = self.dist(sphere2)

for d in range (2):

self.force[d] += (dist[d] * self.gravConst *

(self.mass * sphere2.mass)

/sqrt(dist [0]**2 + dist [1]**2)**3)




UML-Diagram (WH)

Euler

+ integrate


+ integrate





+ calcForce


AstroInter

+ initVis

n

1

1




Schnittstellenklasse zum Steuern der Simulation

• Eine Klasse, die die bisherigen Teile zusammenfuhrt(SimulationInter)

• Dinge, die getan werden mussen:• Initialisierung (Erstellung von Kugeln)• Erstellung des Integrators• Berechnung der Krafte zwischen den Partikeln• Aufruf des Integrators• Erstellung der Zeichenflache• Methode, mit der in einer Schleife Zeichnen,

Kraftberechnung, Anwendung der Randbedingungen undder Integrator wiederholt aufgerufen werden

• Die Klasse implementiert noch keine konkrete Simulation• Dazu wird ein Subklassen von SimulationInter erstellt:

AstroInter




Notige Module

from spheres import GravSphere , LJSphere

from spherevis import Spherevis

from Tkinter import Tk, Canvas

from math import sqrt , pi

from random import random

Konstruktor __init__(self, dt)

• Legt Liste fur Kugeln (spheres[])• Speichert Zeitschrittweite dt

• Erzeugt ein Canvas zur VisualisierungMethode setIntegrator

• Speichert den ubergebenen Integrator ab




Methode calculateForces

• Setzt zunachst die Kraft auf alle Kugeln auf Null• Lauft dann uber alle Kugelpaare und berechnet die Kraft

Methode initVis

• Initialisiert die Visualisierung• Wird durch Subklasse implementiert, da nur da die

Simulationskoordinaten bekannt sindMethode simLoop

• Fuhrt die Simulation endlos durch• In jeden Schritt: Visualisierung, Kraftberechnung, Integration und

Anwendung der Randbedingungen




Implementierung der Klasse SimulationInter

class SimulationInter(object ):

def __init__(self , dt):

self.spheres = []

self.dt = dt

root = Tk()

self.canvas = Canvas(root , width =800, height =800)

self.canvas.pack()

def setIntegrator(self , integrator ):

self.integrator = integrator

def calculateForces(self):

for sphere1 in self.spheres:

sphere1.resetForce ()



if sphere1 != sphere2:

sphere1.calcForce(sphere2)




Implementierung der Klasse SimulationInter (cont.)

def initVis(self):

print "Subclasses have to implement initVis"

exit (1)

def simLoop(self):

while True:

self.vis.drawSpheres(self.spheres)

for i in range (10):

self.calculateForces ()

self.integrator.integrate(self.spheres)




UML-Diagram (WH)

Euler

+ integrate


+ integrate





+ calcForce


AstroInter

+ initVis

n

1

1




Subklasse AstroInter


• Konstruktor der Vaterklasse aufrufen (Zeitschrittweite 1/10 Tag)• Planeten erzeugen

Methode initVis

• Große des sichtbaren Raumausschnitts festlegen• Instanz der Klasse Spherevis erzeugen

class AstroInter(SimulationInter ):

def __init__(self):

SimulationInter.__init__(self , 8640.)

c=3.

self.spheres.append(GravSphere (( 0.00000 , 0.),

(0., 0.), c*20.e9, 1.9891e30 , "yellow")) # Sonne

...

def initVis(self):

self.vis = Spherevis(self.canvas , (-2.e12 ,-2.e12),

( 2.e12 , 2.e12))




Restliche Planeten

self.spheres.append(GravSphere (( 57.9e9, 0.),

(0. ,47872.) , c*10.e9, 3.302 e23)) # Merkur


(0. ,35021.) , c*10.e9, 4.8685 e24)) # Venus


(0. ,29786.) , c*10.e9, 5.9736e24 , "blue")) # Erde


(0. ,24131.) , c*10.e9, 6.4185 e23)) # Mars


(0. ,13070.) , c*10.e9, 1.899 e27)) # Jupiter

self.spheres.append(GravSphere ((1.429e12 , 0.),

(0., 9672.) , c*10.e9, 5.6846 e26)) # Saturn


(0., 6835.) , c*10.e9, 8.6832 e25)) # Uranus


(0., 5478.) , c*10.e9, 1.0243 e26)) # Neptun




Starten der Simulation

• Alle notigen Klassen wurden implementiert• Nun mussen nur noch die konkreten Instanzen erzeugt werden

>>> # notige Module einbinden

>>> from integrator import Euler

>>> from simulationinterface import AstroInter

>>>

>>> # Interface fur Astro -Simulation erstellen

>>> simInter = AstroInter ()

>>> # Integrator mit Zeitschrittweite fur Astro -Sim.

>>> simInter.setIntegrator(Euler(simInter.dt))

>>> # Init Vis. (Abbildung Sim -Koord. auf Canvas)

>>> simInter.initVis ()

>>> # Simulation starten

>>> simInter.simLoop ()




Screenshot der Astro-Simulation




Mogliche Erweiterung?

Welche anderen physikalischen Phanomene lassen sich ahnlichbeschreiben?




Simulation von Molekulen

• Phanomen: Molekule stoßen sich ab bzw. ziehen sich an(klassische Molekulardynamik)

• Basis:• Newtonsche Bewegungsgleichungen:

Fi = mi · ai = mi · ri• Molekule als Partikel betrachtet:

Fi =∑

j Fij• Interaktion zwischen Molekulen: potentielle Energie:

Fij = −∇U(rij )

• Modell: Lennard-Jones-Potential (Kombination ausPauli-Repulsion und van-der-Waals-Anziehung)




Lennard-Jones-Potential

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

0 0.5 1 1.5 2 2.5 3

pote

ntia

l U

distance r

soft sphere potentials

soft sphereLennard−Jonesvan der Waals

σ

ULJ(rij)

=

4ε((

σrij

)12−(σrij

)6)




Simulation von Molekulen III

Anmerkungen• Numerik: analog zur Astro-Simulation• Simulationsprogramm

• Im wesentlichen wie bei der Astro-Simulation• Wir erweitern das bisherige Programm so, dass damit

sowohl Planeten als auch Molekule simuliert werden konnen• Zusatzlich muss das Gebiet begrenzt werden

• Visualisierung wieder mit Tkinter

• Kosten eines Zeitschritts: O(N ∗ N), daher nur wenige Molekulemoglich

• Aufgrund der kurzreichweitigen Kraft (LJ-Potential fallt schnellab) sollten nur Nachbarn von Molekulen zur Kraftberechnungverwendet werden (O(N))




Erweiterung auf MD-Simulation: UML

Euler

+ integrate


+ integrate



Boundary- dim, pos, side

- distance- particleTouches+ boundaryAction

SimulationInter+ spheres+ boundaries+ canvas, dt

- calculateForces- applyBoundaries+ setIntegrator+ initVis+ simLoop

GravSphere

+ calcForce


LJSphere

+ calcForce

AstroInter

+ initVis

MDInter

+ initVis

n

n 1

1




Subklasse LJSphere der Klasse Sphere

• Methode calcForce wird uberschrieben, der Rest von derBasisklass ubernommen

• Lennard-Jones-Potential: U(rij ) = 4ε((σr )12 − (σr )6)

• Zur Vereinfachung: ε = σ = 1• Kraft 1D: F (rij ) = −24(2 · ( 1

r )13 − ( 1r )7)

• Kraft nD: ~F (~rij ) = −24(2 · ( 1d )14 − ( 1

d )8) ·~rij

class LJSphere(Sphere ):


dist = self.dist(sphere2)

absdist = sqrt(dist [0]**2 + dist [1]**2)

for d in range (2):

self.force[d] += 24*(2* absdist **-14

- absdist **-8) * dist[d]




Randbedingungen

• Bisher ist das Simulationsgebiet prinzipiell unendlich groß• Die Molekule wurden in alle Richtungen davon fliegen• Deswegen wird das Gebiet begrenzt• das Verhalten am Rand des Gebiets muss definiert werden:

• Aus- bzw. einstromende Rander• Periodische Rander (was links rausfliegt kommt rechts

wieder rein• Reflektierende Rander (Einfach, deswegen nehmen wir

diese)• Heizende/Kuhlende Rander• ...

Reflektierende Rander• Sobald ein Partikel den Rand beruhrt, prallt es ab• Die Geschwindigkeit senkrecht zur Wand wird invertiert• Die Strecke, die schon uberlappt, wird neuer Abstand zur Wand




Erweiterung auf MD-Simulation: UML (WH)

Euler

+ integrate


+ integrate







GravSphere

+ calcForce


LJSphere

+ calcForce

AstroInter

+ initVis

MDInter

+ initVis

n

n 1

1




Klasse Boundary


• Parameter dim: gibt die Raumrichtung an (0: x, 1: y)• Parameter pos: gibt die Position der Ebene an• Parameter side: -1 fur linke/untere Seite, +1 fur rechts/oben

Methode distance

• Parameter p: Partikel• Berechnet den Abstand des Partikels zum Rand

Methode particleTouches

• Parameter p: Partikel• Gibt True zuruck, wenn das Partikel den Rand beruhrt (Abstand

kleiner als Radius)Methode boundaryAction

• Parameter p: Partikel• Fuhrt bei Beruhrung die Auswirkung des Aufpralls aus




Implementierung der Klasse Boundary

class Boundary(object ):

def __init__(self , dim , pos , side):

self.dim = dim # 0 oder 1

self.pos = pos # beliebiger float -Wert

self.side = side # -1 oder 1

def distance(self , p):

return p.pos[self.dim] - self.pos

def particleTouches(self , p):

return abs(self.distance(p)) <= p.rad

def boundaryAction(self , p):

if self.particleTouches(p):

p.pos[self.dim] -= 2*( self.distance(p)

+p.rad*self.side)

p.vel[self.dim] *= -1




Erweiterung auf MD-Simulation: UML (WH)

Euler

+ integrate


+ integrate







GravSphere

+ calcForce


LJSphere

+ calcForce

AstroInter

+ initVis

MDInter

+ initVis

n

n 1

1




Subklasse MDInter


• Parameter: Anzahl Partikel und Dichte• Konstruktor der Vaterklasse aufrufen (Zeitschrittweite 0.02)• Molekule erzeugen

• Zufallige Positionen schlecht, da sich Molekule dann evtl.uberlappen (hohe Krafte und damit numerische Probleme)⇒ Positionen auf einem Gitter

• Anzahl Partikel so runden, dass ein quadratisches Gittererzeugt werden kann

• Geschwindigkeiten prop. zur Temperatur. Wir versuchennicht, genau zu sein, und wahlen deswegen willkurlichgleichverteilte Geschwindigkeiten zwischen −0.25 und 0.25

• Randbedingungen erzeugen (refl. Rander auf allen Seiten)Methode initVis

• Raum um das Molekulgitter als sichtbaren Bereich festlegen




class MDInter(SimulationInter ):

def __init__(self , numParticles , density ):

SimulationInter.__init__(self , 0.02)

numPerDim = round(sqrt(numParticles ))

self.length = sqrt(pi*numPerDim **2/ density)

for ix in range(int(numPerDim )):

for iy in range(int(numPerDim )):

x = (ix +0.5)/ numPerDim*self.length

y = (iy +0.5)/ numPerDim*self.length

vx = (random () -0.5)/2.

vy = (random () -0.5)/2.

self.spheres.append(LJSphere ((x, y), (vx, vy),

1.0, 1.0, None))

self.boundaries.append(Boundary(0, 0.0, -1))

self.boundaries.append(Boundary(0, self.length , 1))

self.boundaries.append(Boundary(1, 0.0, -1))

self.boundaries.append(Boundary(1, self.length , 1))

def initVis(self):

self.vis = Spherevis(self.canvas , (0.0, 0.0),

(self.length , self.length ))




Evaluation!




Wenn wir eine Simulation auswahlen wollen?

Euler

+ integrate


+ integrate







GravSphere

+ calcForce


LJSphere

+ calcForce

AstroInter

+ initVis

MDInter

+ initVis

n

n 1

1




Einbinden der neuen Klassennotige Module

import argparse

from integrator import Euler

from simulationinterface import AstroInter , MDInter

Option-Parser

# Optionen festlegen

parser = argparse.ArgumentParser(description=’Astro and MD simulation ’)

parser.add_argument(’-n’, ’--numParts ’, dest = "numParts", type=int ,

help=’Number of particles ’, default = 30)

parser.add_argument(’-t’, ’--simType ’, dest=’simType ’, action=’store ’,

help="’Astro ’ or ’MD’ simulation")

parser.add_argument(’--rho’, dest=’rho’, action=’store ’, type=float ,

help="Density", default = 0.1)

options = parser.parse_args ()




Simulationsinterface erstellen

# Simlationsinterface erstellen

if options.simType == None:

print "Bitte Simulationstyp angeben"

exit (1)

if options.simType == "Astro":

simInter = AstroInter ()

elif options.simType == "MD":

n = options.numParts

rho = options.rho

simInter = MDInter(n, rho)




Screenshot der MD-Simulation




Teil XI

Datenstrukturen: Baume, Stacksund Queues




Stacks (Kellerspeicher/Stapel)• Funktioniert wie ein naturlicher Stapel

(z.B. Papierstapel auf dem Schreibtisch)• Elemente werden immer oben auf den

Stapel gelegt (PUSH)• Es kann immer nur das oberste Element

entfernt werden (POP)• Funktionsweise: LIFO (Last In, First Out)

Stacks in Python• Es gibt keine eigene

Stack-Datenstruktur in Python• Listen konnen aber als (effiziente)

Stacks verwendet werden• PUSH: list.append(x)• POP: x = list.pop()

push pop

Sta

ckp

oint

er

pop

pu

sh




Stack-Beispiel: Turme von Hanoi

A B Hilfe

def hanoi(n, a, b, hilfe):

if(n > 0):

hanoi(n-1, a, hilfe , b)

b.append(a.pop())

printhanoi ()

hanoi(n-1, hilfe , b, a)

Randoffscher Algorithmus:

• rekursiv• n=# zu versch. Scheiben, a=woher, b=Zwischenziel, hilfe=wohin




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H

hanoi(1,B,H,A)

B -> H




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H

hanoi(1,B,H,A)

B -> H

A -> B




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H

hanoi(1,B,H,A)

B -> H

A -> B

hanoi(2,H,B,A)

hanoi(1,H,A,B)

H -> A




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H

hanoi(1,B,H,A)

B -> H

A -> B

hanoi(2,H,B,A)

hanoi(1,H,A,B)

H -> A

H -> B




A B Hilfe


if(n > 0):


b.append(a.pop())


hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H

hanoi(1,B,H,A)

B -> H

A -> B

hanoi(2,H,B,A)

hanoi(1,H,A,B)

H -> A

H -> B

hanoi(1,A,B,H)

A -> B




A B Hilfe


if(n > 0):


b.append(a.pop())


Optimale Zugfolge: 2n − 1 Zuge⇒ 34 Jahre bei n=30 und 1

Scheibe/s⇒ 585 Mrd. Jahre bei n=64

hanoi(3,A,B,H)

hanoi(2,A,H,B)

hanoi(1,A,B,H)

A -> B

A -> H

hanoi(1,B,H,A)

B -> H

A -> B

hanoi(2,H,B,A)

hanoi(1,H,A,B)

H -> A

H -> B

hanoi(1,A,B,H)

A -> B




Queues/Warteschlangen

• Funktioniert wie eine naturlicher Warteschlange (z.B.Supermarkt, Postamt, Druckauftrage)

• Elemente werden immer ans Ende der Queue angefugt (PUSH)• Es wird immer nur vom Anfang der Queue entfernt (POP)• Funktionsweise: FIFO (First In, First Out)

Queues in Python• Listen konnen nicht als Queues verwendet werden, da der

Aufwand fur das Entfernen des ersten Elements O(N) ist• Das Module Queue stellt Queues zur Verfugung

>>> from Queue import Queue

>>> q = Queue()

>>> q.put (1); q.put (2); q.put (3)

>>> q.get()

1




Graphen• Wir ersparen uns eine formale Definition!• Ein Graph besteht aus einer Menge Knoten V• Zwischen den Knoten gibt es Kanten E ⊆ V × V

Einsatzgebiete• Verkehrs- und sonstige Netze (kurzeste Wege)• Ablaufplane• Partitionierungsalgorithmen• Viele algorithmische Probleme lassen sich auf Graphen

zuruckfuhren• ...

Was prinzipiell moglich ist• Kanten konnen gerichtet oder ungerichtet sein• Kanten und Knoten konnen Gewichte haben• zwischen zwei Knoten konnen beliebig viele Kanten sein• Ein Graph kann zyklisch, planar oder zusammenhangend (und

jeweils auch das Gegenteil sein)• u.v.m.




Beispiel: Routensuche

Abbildung: Dijkstra bzw. A*; Quelle: Bungartz, Zimmer, Buchholz, Pfluger:Modellbildung und Simulation - Eine anwendungsorientierte Einfuhrung




Datenstrukturen fur GraphenAdjazenzmatrix• Matrix A mit Dimension |V | × |V |• A[i , j] = 1 wenn (i , j) ∈ E• A[i , j] = 0 sonst• Gut geeignet fur dichte Graphen

Adjazenzliste• Liste A mit |V | Listen• n-te Liste aus A enthalt alle Knoten, die uber eine Kante mit

Knoten n verbunden sind• Gut geeignet fur dunne Graphen• In Python alternativ als dictionary realisierbar

Objektorientiert• Knoten sind Objekte• Jeder Knoten enthalt eine Liste mit benachbarten Knoten




(gerichtete) Baume• Gerichteter, zusammenhangender, planarer und azyklischer

Graph• Hat einen Wurzelknoten (nur ausgehended Kanten)• Wenn eine Kante von A zu B gerichtet ist, nennt sich A

Vaterknoten von B bzw. B Kindknoten von A• Knoten ohne Kinder heißen Blatter• Knoten konnen so in einer Hierarchie angeordnet werden, dass

jeder Knoten (außer dem Wurzelknoten) genau eineEingangskante hat, die von einem Knoten in derdaruberliegenden Hierarchieebene liegt.

Einsatzgebiet von Baumen• Speicherung (sortierter) Daten (logarithmischer Zugriff)• Entscheidungs- und Suchbaume• (hierarchische) Gebietszerlegungen• Viele algorithmische Probleme lassen sich auf Baume

zuruckfuhrenT. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Beispiel Breiten- und Tiefensuche

1

2 3

4 75 6

8 9

1

7 2

9 38 4

6 5

Breitensuche Tiefensuche




Tiefensuche in (unsortierten) Baumen/Graphen• Erster Kind-/Nachbarknoten, dann dessen Kind/Nachbar, ...• Worst-Case Laufzeit O(|V |+ |E |)• Es ist moglich, dass der Knoten nicht gefunden wird (Bei

unendlichen Graphen oder wenn Zyklen nicht uberpruft werden)Iterativer Algorithmus mit Stapel

def tiefensuche(startknoten , zielknoten ):

startknoten.besucht = True

stapel = [startknoten]

while len(stapel )>0:

knoten = stapel.pop()

if knoten == zielknoten:

return knoten

for nachbar in knoten.nachbarn:

if nachbar.besucht == False:

nachbar.besucht = True

stapel.append(nachbar)




Breitensuche in (unsortierten) Baumen/Graphen• Erst Knoten mit Entfernung 1 zum Startknoten, dann 2, ...• Bei Baumen: ebenenweiser Durchlauf• Worst-Case Laufzeit O(|V |+ |E |)• Wenn der gesuchte Knoten existiert, wird er auch gefunden

Iterativer Algorithmus mit Warteschlange

def breitensuche(startknoten , zielknoten ):

startknoten.besucht = True

warteschlange = deque([ startknoten ])

while len(warteschlange )>0:

knoten = warteschlange.popleft ()

if knoten == zielknoten:

return knoten

for nachbar in knoten.nachbarn:

if nachbar.besucht == False:

nachbar.besucht = True

warteschlange.append(nachbar)




Suchbaume

• Die bisherigen Verfahren sind nicht effizient (O(|V |+ |E |))• Beide Verfahren ”uninformiert“, d.h. alles wird blind abgesucht• Heuristiken konnen Suche beschleunigen (z.B. A*-Algorithmus

bei Suche nach kurzesten Wegen)• Unterliegen die Elemente des Baumes einer Totalordnung, geht

es mit Suchbaumen noch viel schneller• Suchbaume unterstutzen drei Operationen:

• Einfugen eines Elements (je nach Baum evtl. O(1))• Loschen eines Elements (je nach Baum evtl. O(1))• Suchen eines Elements (O(Baumhoehe))

• Unbalancierte Suchbaume: Baumhohe max. N• Balancierte Suchbaume (z.B. AVL-Baum): Baumhohe = log(N)




Binarer Suchbaum• Erstes hinzugefugtes Element bildet

Wurzel;• Neue Elemente werden als Blatter in

den Baum aufgenommen.• Sie werden von der Wurzel gemaß dem

Sortierkriterium bis zur entsprechendenStelle durchgereicht.

• Die Baumstruktur wird beim Hinzufugenvon Elementen nicht verandert.

• Suche lauft analog zum Einfugen(Vergleich des gesuchten Elements mitdem aktuellen Knoten, dann evtl.Abstieg in linken/rechten Teil, ...);

• Loschen evtl. etwas komplexer(Teilbaume werden umgehangt).

10

5 14

2 8

optimaler Fall

Worst Case

3

12 17

13

10

5

14

2

8

3

12

17

13




AVL-Baum• Binarer Suchbaum, bei dem fur jeden Knoten k gilt, dass die

Hohe von linkem (Hl (k) und rechtem (Hr (k)) Teilbaum sich ummaximal 1 unterscheiden.

• Hohe im schlimmsten Fall ca. 1.44 · log2(N + 2), d.h. Sucheneines Elements immer in O(log(N));

• Fur jeden Knoten gibt es ein Balance-Kriteriumbal(k) = Hr (k)− Hl (k), das immer ∈ −1,0,1 sein muss.

10

5 14

2

normaler Fall

3

12

11

8 17

15 21

16

97

1

0

1

0

1

-1 -1

0

0

10 00

0

14

12

8 17

15 21

7 10

5

2

3

11

-1

-1

-1-1

-1 -1

-1 00

0

0

0

Worst Case




Einfugen eines Elements• Vor dem Einfugen sind alle Balancewerte ∈ −1,0,1.• Zunachst wir der Einfugepunkt gesucht.

Einfugen bei Knoten mit einem Kind• Der bisherige Balancewert des Knotens ist ∈ −1,1.• Der neue Balancewert ist 0.• ⇒ keine Rebalancierung notwendig

Einfugen bei Blattknoten• Der bisherige Balancewert des Knotens ist 0.• Neuer Balancewert ist ∈ −1,1 (aktueller Teilbaum hoher).• Information muss rekursiv nach oben weitergegeben werden.• Bei jedem Vaterknoten wird der neue Balancewert berechnet.

• Bei Balancewert ∈ −1,1 geht man weiter nach oben.• Bei Balancewert 0 kann der Aufstieg abgebrochen werden.• Bei Balancewert ∈ −2,2 muss neu balanciert werden.• Nach Balancierung muss nicht weiter nach oben gegangen

werden, die ubrigen Balancewerte andern sich nicht.T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Balancieren des Baums: Einfachrotationx

y

2

1(bzw. 0)

x

y0

(bzw. -1)

0(bzw. 1)

• Rechter Sohn von x ist y• Linker Teilbaum von x ist am niedrigsten• Linker Teilbaum von y ist nicht hoher als rechter Teilbaum von y• ⇒ y wird Wurzel, x linker Sohn, bisheriger linker Teilbaum von y

wird rechter Teilbaum von x, Rest bleibt gleich• Hohe des neuen Baums gleich wie vor dem Einfugen,

Balancierung beendet• Analoges Vorgehen im spiegelverkehrten Fall




Balancieren des Baums: Doppelrotationx

y2

-1

x2

2/1z

z y

0

0/1

z

yx

1/0/-10/1

-1/0

• Nun ist der linke Teilbaum von y (mit Wurzel z) hoher als derrechte.

• ⇒ Einfachrotation wurde Balance von x von 2 auf -2 andern• Zunachst Rotation, um dafur zu sorgen, dass der rechte

Teilbaum der hochste ist;• Dann eine weitere Einfachrotation wie auf der vorigen Folie.




Loschen eines Elements• Vor dem Loschen sind alle Balancewerte ∈ −1,0,1.• Zunachst wir das zu loschende Element k gesucht

Knoten k ist ein Blatt• Der Knoten wird geloscht, der Teilbaum um 1 niedriger.• ⇒ Balance des Vaters (rekursiv) anpassen• ⇒ evtl. Rebalancierung

Knoten k hat genau einen Sohn s• s tritt an Stelle von k, Teilbaum wird um 1 niedriger.• ⇒ rekursiv Balance anpassen

Knoten k hat zwei Sohne• Bei Balancewert ∈ −1,1 geht man weiter nach oben.• Bei Balancewert 0 kann der Aufstieg abgebrochen werden.• Bei Balancewert ∈ −2,2 muss neu balanciert werden.• Beim Loschen kann sich im Gegensatz zum Einfugen die

Balance des Vaterknotens eines gerade balancierten Teilbaumsauch andern, es muss weiter nach oben gegangen werden.




Baumbasierte KI fur TicTacToeSchnittstelle• Bei jeden Zug wird eine Prozedur aufgerufen, bei der der

KI-Spieler ein ’o’ setzen muss• KI-Spieler erhalt Liste, die das 3x3-Gitter (row-order) enthalt• Drei Zeichen moglich:

• ’ ’: nicht belegt• ’x’: vom Mensch belegt• ’o’: vom PC belegt

• KI muss ein leeres Element (’ ’) mit ’o’ belegen, z.B.:• Vorher: [’o’, ’x’, ’ ’, ’ ’, ’x’, ’ ’, ’ ’, ’ ’, ’ ’]

• Nachher: [’o’, ’x’, ’ ’, ’ ’, ’x’, ’ ’, ’ ’, ’o’, ’ ’]

• Keine Uberprufungen (Bin ich uberhaupt dran, ...)

def pcSchritt(zellen ):

...

zellen [...] = ’o’




Vorgehensweise• Aufbau eines kompletten Entscheidungsbaums (rekursiv)• Jeder Knoten entspricht einer moglichen Spielsituation (z.B.

[’o’, ’x’, ’ ’, ’ ’, ’x’, ’ ’, ’ ’, ’o’, ’ ’])• Die Kinder des Knotens sind die moglichen Zuge des PC.• Fur jeden Knoten bestimmen, welcher Spieler im kompletten

Teilbaum wie oft gewinnt;• Teilbaum mit besten Gewinnchancen auswahlen und das

entsprechende Feld belegen;• Achtung: Vorgehensweise ist weder effizient noch optimal!




Aufbau des Baumes: Parameter fur Konstruktor• self

• zellen: Liste mit der Spielsituation• winproc: Prozedur, die fur die Spielsituation den Gewinner ∈’ ’, ’x’, ’o’ ermittelt

• player=’o’: Aktueller Spieler (fur Wurzelknoten immer ’o’)• zug=None: Zug der zur aktuellen Spielsituation gefuhrt hat (fur

Wurzelknoten None)Aufbau des Baumes: member-Variablen• zellen: Liste mit der Spielsituation• kinder: Liste der Kinder des Baumes• zug: Welcher Zug (0-8) hat zu der Spielsituation gefuhrt• gewinner Evtl. gibt es fur die Situation schon einen Gewinner (’x’

oder ’o’), sonst ’ ’




Konstruktor der Klasse

class TicTacToeBaum(object ):

def __init__(self , zellen , winproc , player=’o’,

zug =0):

self.zellen = zellen

self.kinder = []

self.zug = zug

nextplayer = ’x’: ’o’, ’o’: ’x’

self.gewinner = winproc(zellen)

if self.gewinner == ’ ’:

for i in range (9):

if (zellen[i] == ’ ’):

copy = list(zellen)

copy[i] = player

self.kinder.append(TicTacToeBaum(copy ,

winproc , nextplayer[player], i))




Gewinner fur Teilbaume bestimmen• Erneut rekursive Berechnung• dictionary anzGewinne = ’ ’: 0, ’x’: 0, ’o’: 0

• Wenn ein Knoten keine Kinder hat, ist anzGewinne[gewinner] = 1

• Sonst wird die Summe der Gewinne der Kinder ermittelt

def zaehleGewinner(self):

self.anzGewinne = ’ ’: 0, ’x’: 0, ’o’: 0

if len(self.kinder )==0:

self.anzGewinne[self.gewinner] += 1

for kind in self.kinder:

(nowin , xwin , owin) = kind.zaehleGewinner ()

self.anzGewinne[’ ’] += nowin

self.anzGewinne[’x’] += xwin

self.anzGewinne[’o’] += owin

return (self.anzGewinne[’ ’], self.anzGewinne[’x’],

self.anzGewinne[’o’])




Wahl des ”optimalen“Feldes• Einschrankung:

• Perfekte Strategie relativ aufwandig• Hier deswegen nur eine mittelmaßige (damit es auf eine

Folie geht)• Fur alle Spielzuge (alle Kinder) wird die ”Gewinnchance“

bestimmt.• Wenn in einer Spielsituation nur noch der PC-Spieler gewinnen

kann (auch kein Unentschieden), ist das optimal, dieGewinnchance riesig.

• Wenn zusatzlich noch Unentschieden moglich sind, ist dieChance immer noch groß.

• Wenn beide gewinnen konnen, wird die Chance als Quotient derbeiden Haufigkeiten gewahlt.

• Wenn nur noch der Mensch gewinnen kann, entspricht dieChance der negativen Haufigkeit der menschlichenGewinnmoglichkeiten.

• Das Kind (der Zug ∈ 0− 8) mit der hochsten Chance (nachobiger Heuristik) wird ausgewahlt.




def optWahl(self):

chance = -1e100

wahl = self.kinder [0]. zug

for kind in self.kinder:

nowin = float(kind.anzGewinne[’ ’])

xwin = float(kind.anzGewinne[’x’])

owin = float(kind.anzGewinne[’o’])

if owin > 0 and xwin == 0 and nowin == 0:

neueChance = 1e100

elif owin > 0 and xwin == 0:

neueChance = 1e10*owin/nowin

elif owin > 0 and xwin > 0:

neueChance = owin/xwin

else:

neueChance = -xwin

if neueChance > chance:

chance = neueChance

wahl = kind.zug

return wahl




Todo: Summary-Folie als Rausschmeisser!




Teil XII

Wissenschaftliches Rechnen inPython




Nochmal Modul math

• Konstanten pi und e

• Funktionen fur int und float

• Alle Ruckgabewerte sind float

ceil(x)

floor(x)

exp(x)

fabs(x) # wie abs(), nur immer float

ldexp(x, i) # x * 2**i

log(x [,base])

log10(x) # == log(x, 10)

modf(x) # (Nachkommateil , Integerteil)

pow(x, y) # x**y

sqrt(x)




Trigonometrische Funktionen (Bogenmaß)

cos(x); cosh(x); acos(x)

sin(x); ...

tan(x); ...

degrees(x) # rad -> deg

radians(x) # deg -> rad

inf/nan

float("inf")

float("-inf")

float("nan")

isinf(x) # Uberprufung auf Unendlichkeit

isnan(x) # Uberprufung auf NaN

Complexe Zahlen: cmath




Und RICHTIGE Mathematik?

Bisherige Sequenztypen (list, tuple, . . . )• Konnen als Arrays verwendet werden• Konnen verschiedene beliebige Typen enthalten

• Sehr flexibel, aber auch langsam• Schleifen sind nicht sehr effizient

• Vektoren/Matrizen und Operationen auf diesen sind extremumstandlich

• Fur effizientes wissenschaftliches Rechnen sind andereDatentypen und Methoden notig

Module• NumPy• Matplotlib• SciPy




Modul numpy

Homogene Arrays• NumPy stellt beliebig-dimensionale Arrays bereit• Alle Elemente des Arrays haben den gleichen Typ• Beispiel

from numpy import *

a = array ([[1 ,2 ,3] ,[4 ,5 ,6]])

print a

type(a)

a.shape

print a[0,2]

a[0,2] = -1

b = a*2

print b




Homogene Arrays (2)• Arrays konnen aus (verschachtelten) Sequenztypen erstellt

werden• Direkter Zugriff mit []

a = array ([1,2,3,4,5,6,7,8])

a[1]

a = array ([[1,2,3,4],[5,6,7,8]])

a[1,1]

a = array ([[[1 ,2] ,[3 ,4]] ,[[5 ,6] ,[7 ,8]]])

a[1,1,1]

• Eigenschaften von Arrays

a.ndim # Anzahl Dimensionen

a.shape # Anzahl Eintrage in jeder Dimension

a.size # Anzahl Elemente

a.dtype # Datentyp der Elemente

a.itemsize # Notige Anzahl Bytes fur den Datentyp




Datentypen

• Exakte Typen (ahnlich C/C++) konnen angegeben werden• Ohne Angabe werden Standardtypen verwendet• Einige Typen (v.a. unterschiedlicher Speicherbedarf):

• int, int8, int16, int32, int64,• float, (float8, float16,) float32, float64,• complex, complex64,• bool, character, object

array ([[1,2,3],[4,5,6]], dtype=int)

array ([[1,2,3],[4,5,6]], dtype=complex)

array ([[1,2,3],[4,5,6]], dtype=int8)

array ([[1 ,2 ,3] ,[4 ,5 ,1000]] , dtype=int8) # falsch

array ([[1,2,3],[4,5,"hi"]], dtype=object)




Numerische Typen - ganze Zahlen (int)

• Darstellung im Rechner: zyklisch (neg,0,pos)• Beispiel overflow (C++, python):

#include <iostream >

#include <limits >

int main()

int d=std:: numeric_limits <int >:: max ()-10;

while(d > 0) std::cout << d++ << std::endl;

return 1;

import sys

d = sys.maxint -10

while d>0:

d += 1

print d




Arrays initialisieren• Standard-Matrizen (optionaler Parameter: dtype)

arange ([a,] b [,stride ]) # analog zu range , 1D

zeros( (3,4) )

ones( (1,3,4) )

empty( (3,4) ) # keine Init. (schnell)

linspace(a, b [, n]) # n aquid. Werte in [a,b]

logspace(a, b [, n]) # log zw. 10^a und 10^b

identity(n) # Identitat (Diagonalmat .)

fromfunction(lambda i,j: i+j, (3,4), dtype=int)

def f(i,j):

return i+j

fromfunction(f, (3,4), dtype=int)




Arrays manipulieren

>>> a = arange (12); a

array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])

>>> b = a.reshape ((3 ,4)); b; a

array ([[ 0, 1, 2, 3],

[ 4, 5, 6, 7],

[ 8, 9, 10, 11]])

array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])

>>> a.resize ((3 ,4)) # a wird verandert

>>> a.transpose () # a wird nicht verandert

>>> a.flatten () # a wird nicht verandert

# Beispiel:

a = arange (144)

a.resize ((12 ,12))




Noch mehr Manipulations- und Initialisierungsmethoden• Neue Matrizen aus existierenden Matrizen erstellen

a = arange (12); a.resize ((3 ,4))

copy(a)

diag(a); tril(a); triu(a)

empty_like(a) # Kopiert nur die Form

zeros_like(a)

ones_like(a)

a = loadtxt("matrix.txt") # fromfile () bei Binardatei

• Matrizen ausgeben

a.tolist ()

savetxt("matrix.txt", a) # tofile () bei Binardatei




Lesen und Schreiben von Array-Elementen

• Die von Listen bekannte Indizierung kann verwendet werden• Dimensionen werden mit Komma separiert


a[1]

a[1:3 ,:]

a[: ,::2]

a[::2 ,::2]

a[::2 ,::2] = [[0, -2, -4],[-10, -12, -14]]

a[1::2 ,1::2] = -1*a[1::2 ,1::2]

• Selektiver Zugriff

a[a > 9]

a[a > 9] *= -1




Uber Elemente von Matrizen iterieren


for row in a:

print row

b = arange (30); b.resize ((2,3,4))

for row in b:

for col in row:

print col

for entry in a.flat: # Iterator

print entry




Rechnen mit Arrays• Schnelle built-in Methoden fur Arrays


3*a

a**2

a+a**2

sin(a)

sqrt(a)

prod(a)

sum(a)

it = transpose(a)

x = array ([1,2,3])

y = array ([10 ,20 ,30])

inner(x, y)

dot(it , x)

cross(x,y)




Rechnen mit Arrays

• Es gibt noch viel mehr. . .

var() cov() std()

mean() median ()

min() max()

tensordot ()

...

• Matrizen (mit mat) sind Unterklassen von ndarray, aber striktzweidimensional, mit zusatzlichen Attributen

m = mat(a)

m.T # Transponierte

m.A # Als 2d Array

m.H # Konjugiert Transponierte




Untermodule

Modul numpy.random• Zufallszahlen aus vielen unterschiedlichen Verteilungen• Machtiger als das Standard-Modul random• Arbeitet auf Arrays und gibt Arrays zuruck

from numpy.random import *

binomial (10, 0.5) # 10 Versuche , Erfolg 50%

binomial (10, 0.5, 15) # 15 Mal voriges

binomial (10, 0.5, (3 ,4)) # (3x4) Array

randint(0, 10, 15) # [0,10), int

rand() # [0,1), float

rand (3,4) # (3x4) Array




Untermodule (2)

Modul numpy.linalg• Die wichtigsten Werkzeuge fur lineare Algebra

from numpy.linalg import *


norm(a); norm(x) # Matrix - und Vektornorm

inv(a) # Inverse der Matrix

solve(a, b) # LAPACK LU Zerlegung

det(a) # Determinante

eig(a) # Eigenwerte und Eigenvektoren

cholesky(a) # Cholesky -Zerlegung

Modul numpy.fft• Fourier-Transformation

Und weitere . . .




Version Mania

Aktuelle Situation

python

MatplotlibSciPy

NumPy

IPython

pylab




Version Mania

Probleme:• Numpy, scipy, pylab, ipython und matplotlib werden oft

gleichzeitig benutzt• Die Pakete hangen voneinander ab (matplotlib benutzt z.B.

numpy-Arrays• Je nach Betriebssystem(sversion) mussen evtl. unterschiedliche

Pakete installiert werden (D.h. Der Modulname im import-Befehlkann auch unterschiedlich sein).

Vision: Alles in einem (neuen) Modul PyLab vereinigen!• gibt es als inoffizielles Paket• Achtung Name: wieder pylab!




Matplotlib




Matplotlib

Wozu Matplotlib?• Objektorientierte Bibliothek fur zweidimensionale plots• Entworfen mit den Ziel einer ahnlichen Funktionalitat wie Matlab

plots• Insbesondere zum Darstellen wissenschaftlicher Daten, baut auf

numpy-Datenstrukturen auf




Moglichkeiten zum Importierenipython interaktiv

> ipython

import scipy , matplotlib.pylab

x = scipy.randn (10000)

matplotlib.pylab.hist(x, 100)

matplotlib.pylab.show()

> ipython

import numpy.random , matplotlib.pylab

x = numpy.random.randn (10000)

matplotlib.pylab.hist(x, 100)

matplotlib.pylab.show()

ipython mit pylab-Modus

> ipython -pylab

x = randn (10000)

hist(x, 100)T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Beispiel - Erster Plot

http://matplotlib.sourceforge.net/users/screenshots.html

from pylab import *

x = arange (0.0, 2*pi , 0.01)

y = sin(x)

plot(x, y)

xlabel(’Label fuer x-Achse ’)

ylabel(’Label fuer y-Achse ’)

title(’Einfacher sin plot’)

grid(True)

show()




Beispiel - Subplots verwenden

from pylab import *

def f(t):

s1 = cos (2*pi*t)

e1 = exp(-t)

return multiply(s1,e1)

t1 = arange (0.0, 5.0, 0.1)

t2 = arange (0.0, 5.0, 0.02)

t3 = arange (0.0, 2.0, 0.01)

show() # Gibt einen Fehler aber hilft ;-)

subplot (2,1,1) # Zeilen , Spalten , welche anzuzeigen

plot(t1, f(t1), ’go’, t2, f(t2), ’k--’)

subplot (2,1,2)

plot(t3, cos (2*pi*t3), ’r.’)




Beispiel - Subplots verwenden

# Fortsetzung der vorigen Folie

subplot (2,1,1)

grid(True)

title(’A tale of 2 subplots ’)

ylabel(’Gedaempfte Oszillation ’)

subplot (2,1,2)

grid(True)

xlabel(’Zeit (s)’)

ylabel(’Ungedaempft ’)




Beispiel - Histogramm

from numpy import *

from matplotlib.mlab import *

from matplotlib.pyplot import *

#alternativ: ipython -pylab und unten random. weglassen

mu, sigma = 100, 15

x = mu + sigma*random.randn (10000)

n, bins , patches = hist(x, 50, normed=1, \

facecolor=’green’, alpha =0.75)

# Linie mit ’best fit’ einzeichnen

y = normpdf( bins , mu , sigma)

plot(bins , y, ’r--’, linewidth =1)

axis ([40, 160, 0, 0.03])

show()




SciPy




Mehr als NumPy?

• SciPy hangt von NumPy ab• SciPy arbeiten mit NumPy arrays• Stellt Funktionalitat fur Mathematik, Naturwissenschaft und

Ingeneursanwendungen zur Verfugung• Noch in Entwicklung• Bei NumPy geht es im Wesentlichen um (N-dimensionale)

Arrays.• SciPy stellt eine große Zahl an Werkzeugen zur Verfugung, die

diese Arrays verwenden.• SciPy beinhaltet die NumPy Funktionalitat (nur ein import notig).• Es gibt noch einige weitere Bibliothen fur wissenschaftliches

Rechnen (teilweise aufbauend auf NumPy und SciPy).• Hier daher nur eine kurze Ubersicht• Weitere Informationen auf www.scipy.org




SciPy - Unterpaketecluster Clustering Algorithmenconstants Physikalische und mathematische Konstantenfftpack Fast Fourier Transformationintegrate Integration und Loser fur gewohnliche DGLsinterpolate Interpolation und Splinesio Ein- und Ausgaberoutinenlinalg Lineare Algebramaxentropy Maximale Entropiendimage N-dimensionale Bildverarbeitungodr Orthogonal distance regressionoptimize Optimierung und Nullstellensuchesignal Signalverarbeitungsparse Dunne Matrizen und zugehorige Routinenspatial Raumliche Datenstrukturen und Algorithmenspecial Spezielle Funktionenstats Statistische Verteilungen und Funktionenweave Integration von C/C++ Programmen




Spezielle Funktionen

• Airy-Funktion• Elliptische Integralfunktion• Bessel-Funktionen (+ Nullstellen, Integrale, Ableitungen,...)• Struve-Funktion• Jede Menge statistische Funktionen• Gamma-Funktion• Hypergeometrische Funktionen• Orthogonale Polynome (Legendre, Chebyshev, Jacobi,...)• Parabolische Zylinderfunktion• Mathieu-Funktion• Kelvin-Funktion• ...




Beispiel: Interpolation - Linear

import numpy as np

import matplotlib.pyplot as plt

from scipy import interpolate

x = np.arange(0, 2.25*np.pi , np.pi/4)

y = np.sin(x)

f = interpolate.interp1d(x, y)

xnew = np.arange (0 ,2.25*np.pi, np.pi/4)

ynew = f(xnew)

plt.plot(x,y,’o’,xnew ,ynew ,’-’)

plt.title(’Lineare Interpolation ’)

plt.show()




Beispiel: Interpolation - Kubische Splines

import numpy as np


from scipy import interpolate

x = np.arange(0, 2.25*np.pi , np.pi/4)

y = np.sin(x)

spline = interpolate.splrep(x,y,s=0)

xnew = np.arange (0 ,2.02*np.pi,np.pi/50)

ynew = interpolate.splev(xnew , spline)

plt.plot(x,y,’o’,xnew ,ynew)

plt.legend ([’Interp -Punkte ’,’Kubischer Spline ’])

plt.axis ([ -0.05 ,6.33 , -1.05 ,1.05])

plt.title(’Interpolation mit kubischen Splines ’)

plt.show()




Teil XIII

Simulation von PartiellenDifferentialgleichungen (PDE):

Warmeleitungsgleichung




Werbung: 2 Stellen fur Stud. HilfskrafteFerienakademie• Umfang: ca. 5–7h/Woche (ca. EUR 211-295/Monat)• Aufgabe: Weiterentwicklung eines online-Bewerbungssystems• Anforderungen

• Programmierung mit java• Nutzung DB2 Datenbankanbindung• hilfreich, aber kein Muss: UNIX-Basiserfahrung

• Kontakt: Sebastian Rettenberger ([email protected])Serviceburo IN• Umfang: 6-10 h/Woche (ca. EUR 253–422/Monat)• Aufgabe: Weiter- und Neuentwicklung von

MS-Access-Datenbanken• Anforderungen: Einfuhrung in Datenbanken• Umfeld: MS-Access Datenbanken, Windows :-)• Kontakt: Dr. Herbert Ehler ([email protected])




ODE vs. PDE• Differentialgleichungen bei der Molekulardynamik:

• nur eine unabhangige Variable: Zeit• gewohnliche Differentialgleichungen (ODE)

• Bei vielen physikalischen Problemen:• zwei oder mehr unabhangige Variablen: Zeit und bis zu drei

Raumrichtungen• partielle Ableitungen nach allen Variablen, deshalb partielle

Differentialgleichung (partial differential equation, PDE)• Einsatz der Modellwerkzeuge ODE und PDE auch alternativ:

• hohere Genauigkeit bei PDE (raumliche Heterogenitaten)• hoherer Aufwand bei PDE (Diskretisierung von Zeit und

Raum)• Beispiel: Populationsdynamik




Raumauflosende Populationsmodelle

• Einerseits: rein zeitabhangige Modelle manchmal zu grob• Bevolkerungsentwicklung in den USA (1850er, California

gold rush): starke Ost-West-Komponente• Migration der Weltbevolkerung im 21. Jahrhundert• Vorhersage von Heuschreckenplagen in Afrika: flachige

Ausbreitung, diffusive und andere Effekte• deshalb Zielgroße eher p(x , t) oder p(x , y , t) anstelle von p(t)• Erhohung der Komplexitat:

• Modelle: Zusammenhang von Raum und Zeit?• Numerik: Speicher- und Rechenaufwand u.v.m.

• andererseits: hohere Genauigkeit uberhaupt erforderlich?• Europa: viele Auswanderer (egal, wo das Gold ist)• USA: wichtig zu wissen, wo Infrastruktur bereitzustellen




Physikalische Phanomene, die PDEs erfordernStromungsmechanik/Thermodynamik• Wo entsteht ein Tornado?• Ist die gegebene Karosserie aerodynamisch gunstig?

Strukturmechanik• Halt das Gebaude einer Belastung stand?• Wo sind Sollbruchstellen?

Verfahrenstechnik• Wo wird es wie heiß in einem (nuklearen / chemischen /

biologischen) Reaktor?Elektromagnetismus• Wo im Transistor ist die Elektronendichte wie hoch?

Geologie• Wann, wo und wie heftig wird es Erdbeben geben?

...T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Mit oder ohne Zeit?

stationare Phanomene• keine Zeit-, sondern nur Ortsabhangigkeit• bei eindimensionalem Raum: wieder ODE!• Beispiele: Gleichgewichtssituationen

• Auflosung einer Substanz in Wasser ohne außere Einflusse• Temperaturverteilung in konstant beheiztem Raum

instationare Phanomene• Zeit- und Ortsabhangigkeit• auf jeden Fall PDE!• Beispiele:

• Oszillationen: mechanische Belastung einer Schiffschaukel• Turbulenz (vgl. Stromung uber Wehr, Wirbel im Nachlauf

startender Flugzeuge)• Regelkreis (standig angepasste Randbedingungen)




Beispiel: WarmeleitungKernproblem der Thermodynamik• Warme beeinflusse die Oberflache eines Objekts• Ziel: Ausbreitung bzw. Verteilung im Objekt vorhersagen• an der Oberflache: Randbedingungen• Beginn der Simulation: Anfangsbedingungen• Materialeigenschaften (Warmeleitfahigkeit etc.) beeinflussen

WarmeleitungBeispiele• ein Hitzdraht• ein Kochtopf auf einer Herdplatte• Kuhlwasser im Reaktor eines Kernkraftwerks• ein Zimmer im Winter: wo die Heizung platzieren?• ein Zimmer im Sommer: Aufheizung an Fensterflachen

Objekt der Begierde i.W.: Temperatur T

T (x ; t) oder T (x , y ; t) oder T (x , y , z; t)




Szenario: Kochtopf

• Gebiet Ω ⊂ R3 , Rand ∂Ω ,Zeitintervall [0, τ ]

• Zu Beginn: Topf und Herdplatte kalt• Nach dem Einschalten: Platte schnell heiß, Wasser erhitzt sich

langsam; am restlichen Gebietsrand: Kuhlung1D?• Nur die Hohe als Raumdimension (unten heiß, oben kalt)• Aber: Kuhlung am Rand, Platte evtl. kleiner als Topf

2D?• Zusatzlich zur Hohe noch den Abstand zur Hohenachse• Aber: Ist die Herdplatte wirklich rotationssymmetrisch?

3D?• Auflosung des kompletten Raumes• Aber: sehr rechen- und speicheraufwandig




Herleitung des Modells

Die folgenden zwei Modellteile mussen hergeleitet werdenPDE• Beschreibt die Wechselwirkungen von Temperaturanderungen in

Bezug auf Raum und Zeit• Eine Gleichung fur eine gesuchte Funktion ausreichend• Gleichung bestimmt Schar von Losungen

Anfangs- und Randbedingungen• Zur eindeutigen Festlegung der gesuchten Losung• Anfangsbedingungen: Temperaturfeld im gesamten Gebiet zum

Startzeitpunkt (wie bei ODE)• Randbedingungen: Temperaturdaten am Rand des Gebiets zu

allen Zeiten




Instrumentarium: DifferentialoperatorenGradient ∇ (Nabla-Operator)• Angewandt auf skalare Funktion f : Rn → R• Vektor der partiellen Ableitungen• Zeigt in die Richtung des steilsten Anstiegs• Die Lange ist ein Maß fur die Steilheit

Divergenz div• Angewandt auf Vektorfeld ~F : Rn → Rn

• div ~F :=∑n

i=1∂Fi∂xi

• Gibt fur jeden Punkt des Feldes an, ob etwas zu- oder abfließt• Positive Werte: Quellen, Negative Werte: Senken

Laplace-Operator ∆

• Angewandt auf skalare Funktion f : Rn → R• Ist definiert als die Divergenz des Gradienten• ∆f := div∇f =

∑ni=1

∂2f∂x2

i.




Prinzipielle Vorgehensweise

• Start bei einem grundlegenden physikalischen Gesetz:• Der beruhmte Apfel• Die Erhaltung von Masse• Die Erhaltung von Energie oder ...

• Daraus Gewinnung einer Bilanzgleichung (typischerweise einsummatorischer Zusammenhang mit Integralen)

• Zuhilfenahme der Analysis: Vereinfachung• z.B. Gaußscher Integralsatz∫

G div ~F (~x)d~x =∫∂G~F (~x)~n

−→dS

•”Integral verschwindet stets“zieht ”Integrand Null“nach sich

• Vereinfachung durch physikalische Einschrankungen• Keine Quellterme etc.

• Am Ende dann die gewunschte (Differential-) Gleichung




Herleitung der WarmeleitungsgleichungAnnahme: homogenes Material• Dichte ρ konstant• Spezifische Warmekapazitat c konstant (Warmekapazitat

beschreibt, wieviel Energie benotigt wird, eine bestimmte Mengedes Stoffs auf eine bestimmte Temperatur zu bringen)

Warmeanderung• Betrachtung eines endlichen Volumenstucks V• Warme: Q =

∫V cρT (x ; t) dx

• Innerhalb V Energieerhaltung (Sofern keine Quellen undSenken)

• Anderung der Warme nur durch Zu- und Abfluß uber Oberflache• Anderung proportional zur Normalenableitung ∂T

∂n an derOberflache, also proportional zum Skalarprodukt von ∇T und ~n

• Proportionalitatsfaktor: Warmeleitfahigkeit k




Bilanzgleichung• Anderung der Warme entspricht Zu- und Abflussen

dQdt

=

∫

Vcρ∂T∂t

(~x ; t) d~x =

∫

∂Vk∇T (~x ; t)~n

−→dS

• Gaußscher Integralsatz∫

Vcρ∂T∂t

(~x ; t) d~x =

∫

Vk∆T (~x ; t) d~x

• Gilt fur jedes beliebige Volumenstuck

cρ∂T∂t

(~x ; t) = k∆T (~x ; t)

• Mit κ = k/(c%) folgt die Warmeleitungsgleichung:

∂T (~x ; t)∂t

= κ∆T (~x ; t) ,




Partielle DifferentialgleichungenMotivation• Bisher nur Beispiele mit gewohnlichen Differentialgleichungen• Z.B. bei Molekulen: Bewegung jedes Molekuls durch eine ODE

beschrieben, nur Abhangigkeit von der Zeit• Viele physikalische Prozesse lassen sich nur in Abhangigkeit von

mehreren Variablen beschreiben• Haufiger Fall: Zeitabhangigkeit und zwei oder drei

RaumvariablenAllgemeine PDE in zwei Dimensionen

F (x , y ,u(x , y),∂u(x , y)

∂x,∂u(x , y)

∂y,∂2u(x , y)

∂x2 , ...,∂(i+j)

∂x i∂y j ) = 0

Lineare PDE zweiter Ordnung (konstante Koeffizienten)

a∂2u(x , y)

∂x2 + b∂2u(x , y)

∂y2 + ...+ e∂u(x , y)

∂y+ fu(x , y) = g(x , y)




Diskretisierung

• Prinzip ahnlich wie bei ODEs (Euler-Methode)• Behandlung kontinuierlicher Gleichung in endlicher Zeit mit

endlichem Speicher• Zerlegung eines Gebiets (n-dimensional) in eine endliche Zahl

kleinerer Einheiten• Zerlegung wird meist Gitter genannt

Diskretisierung des Raums: Gitter• Je nach Diskretisierungsmethode und Problemstellung

unterschiedliche Gitter• strukturierte Gitter (z.B. kartesisch): Gitterpunkte mussen nicht

explizit gespeichert werden• unstrukturierte Gitter: Koordinaten und Topologie mussen

gespeichert werden• Bei unstrukturierten Gittern ist Adaptivitat leichter umzusetzen




Diskretisierungsansatze fur Gleichungsoperatoren• Finite Differenzen• Finite Volumen• Finite Elemente




Verschiedene Gittertypen

strukturiertes Gitter 2D unstrukturiertes Gitter 2D

adaptives Gitter 2D adaptives Gitter 3D




Adaptivitat

• Feineres Gitter⇒ teurer und genauer• Oft: hohe Genauigkeit nur an wenigen Stellen notig!• Beispiele:

• Strukturmechanik: ”Sollbruchstellen“• Stromungsmechanik: Wirbel• Warmeleitung: dunne Bauteile




Finite Volumen

• Das Gebiet wird in endliche Zahl an Zellen zerlegt• Strukturierte und unstrukturierte Gitter moglich• Fur jede Zelle gilt ein Erhaltungssatz (Warme verschwindet nicht)• Mittlere Temperatur: Integral uber das Volumen der Zelle• Gaußscher Integralsatz: Integral uber den Rand der Zelle• ⇒Warmeanderung nur durch Zu- und Abflusse uber den Rand• Man erhalt fur jede Zelle eine gewohnliche Differentialgleichung• ⇒ Diskretisierung gibt algebraische Gleichung• Fur jede Zelle gibt es daher eine Gleichung• Insgesamt ein lineares Gleichungssystem




Finite Elemente

• Zu kompliziert, um hier en detail erklart zu werden• Gebiet wird in (unstrukturierte) Elemente unterteilt• PDE wird nicht uberall erfullt, sondern nur gemittelt• Letztlich erhalt man fur jeden Freiheitsgrad eine Gleichung• Damit auch wieder ein lineares Gleichungssystem




Finite Differenzen

• Vergleichsweise einfach (strukturierte Gitter, einfache Theorie)• PDE wird direkt diskretisiert: Ersetzung der Ableitungen durch

Differenzenquotienten• Bei erster Ableitung: Approximation durch Sekante• Vorwartsdifferenzenquotient: T ′(x0) = T (x0+h)−T (x0)

h

• Ruckwartsdifferenzenquotient: T ′(x0) = T (x0)−T (x0−h)h

• Beide Male: Differenz zweier Funktionswerte durch Differenz derx-Werte

T

T 0

T 1

h xx0 x1




3-Punkte-Stern• Fur die 1D Warmeleitungsgleichung• Diskretisierung der zweiten Ableitung• Analog zur ersten Ableitung: Differenz zweier Ableitungswerte

durch die Differenz der zugehorigen x-Werte• Die beiden Ableitungswerte erhalt man durch Vorwarts- und

Ruckwartsdifferenzenquotient

T ′′(x0) =T (x0+h)−T (x0)

h − T (x0)−T (x0−h)h

h

=T (x0 + h)− 2T (x0) + T (x0 − h)

h2 .

• Maschenweite h bestimmt wieder Aufwand und Genauigkeit• Veranschaulichung an einem Stab:

0

-2 11

h

T 1 T 2 T i−1 T i T i1 T nT n−1

h L




Anwendung auf die stationare Warmeleitungsgleichung• stationare Warmeleitungsgleichung: κ∆T (x) = 0• Fur jeden der n diskreten Punkte gilt:

κ · Ti+1 − 2Ti + Ti−1

h2 = 0

• Und nach Vereinfachung:

Ti+1 − 2Ti + Ti−1 = 0

• Daraus folgt ein Lineares Gleichungssystem:

−2 1 0 . . . 0

1 −2 1...

0 1 −2. . . 0

.... . . . . . 1

0 . . . 0 1 −2

·

T1T2...

Tn−1Tn

=

−T00...0

−Tn+1

.




Aufbau von LGS fur 1D-Diskretisierung

• Verwende numpy Matrizen• Erzeuge zufallige rechte Seite

import numpy

N = 9

A = numpy.zeros((N, N))

A[0, 0] = -2

A[0, 1] = 1

for i in range(1, N-1):

A[i, i - 1] = A[i, i + 1] = 1

A[i, i] = -2

A[-1, -1] = -2

A[-1, -2] = 1

b = numpy.random.random ((N, 1))

x = numpy.linalg.solve(A, b)




5-Punkte-Stern

• Bei 2D-Gleichung: zwei partielle zweite Ableitungen ∂2T∂x2 und ∂2T

∂y2

• In beiden Dimensionen analog zur 1D-Vorgehensweise

∂2T∂x2 =

T (x + h, y)− 2T (x , y) + T (x − h, y)

h2 ,

∂2T∂y2 =

T (x , y + h)− 2T (x , y) + T (x , y − h)

h2 .

• Anwendung auf κ∆T (x) = 0, Temperatur in Gitterzellen Ti,j

κ

(Ti+1,j + Ti,j+1 − 4Ti,j + Ti−1,j + Ti,j−1

h2

)= 0

• Wieder vereinfachen:

Ti+1,j + Ti,j+1 − 4Ti,j + Ti−1,j + Ti,j−1 = 0 .




Veranschaulichung der 2D-Diskretisierung

-4 1 1110 0 0 0 0 0

j=m

l+1

l

l-1

j=1

i=1 i=nk-1 k k+1




LGS fur 2D-Diskretisierung• Anordnung der Zellen aus dem 2D-Gitter zeilenweise im

1D-Vektor der Unbekannten• In der folgenden Darstellung: 3x3 Gitterzellen• Daher nur eine Zelle nicht am Rand (und damit rechte Seite

gleich Null)

−4 1 11 −4 1 1

1 −4 11 −4 1 1

1 1 −4 1 11 1 −4 1

1 −4 11 1 −4 1

1 1 −4

·

T1,1T1,2T1,3T2,1T2,2T2,3T3,1T3,2T3,3

=

−T0,1 − T1,0−T0,2

−T0,3 − T1,4−T2,0

0−T2,4

−T3,0 − T4,1−T4,2

−T4,0 − T3,4




Demo Comsol (kommerzielles Werkzeug)

Aluminium Heat Sink

Copper Cable

273 K

1e9 W/m^3




Teil XIV

Losung linearerGleichungssysteme




Wo tauchen LGS uberall auf?




Gauß-Algorithmus

• Zwei Schritte: Vorwartselimination und RuckwartssubstitutionVorwartselimination• Erzeugen einer Stufenform• Zeilen durfen mit Konstanten multipliziert werden• Zeilen durfen addiert werden werden

a11 a12 a13a21 a22 a23a31 a32 a33

·

x1x2x3

=

b1b2b3

- 1. Zeile · a21

a11

- 1. Zeile · a31a11

a11 a12 a130 a22 a230 a32 a33

·

x1x2x3

=

b1

b2

b3

- 2. Zeile · a32a22

a11 a12 a130 a22 a230 0 a33

·

x1x2x3

=

b1

b2

b3

- 2. Zeile · a32a22




Ruckwartssubstitution• LGS muss in Stufenform (obere Dreiechsmatrix) vorliegen• Die unterste Gleichung hat nur eine Variable• Deren Wert kann direkt berechnet werden• Zweitunterste Zeile hat zwei Variablen, eine davon unbekannt• . . .• nte Zeile von unten hat n Variablen, eine unbekannt

x3 = b3/a33

x2 = b2 − a23 · x3/a22

x1 = b1 − a12 · x2 − a13 · x3/a11

• Zeile in einer großeren Matrix:

xi = (bi −n∑

j=i+1

aij · xj )/aii




Beispiel

2 3 14 7 32 4 4

·

x1x2x3

=

122

2 3 10 1 10 1 3

·

x1x2x3

=

101

2 3 10 1 10 0 2

·

x1x2x3

=

101

x3 = 1/2 = 0.5x2 = (0− 1 · 0.5)/1 = −0.5x1 = (1− 3 · (−0.5)− 1 · 0.5)/2 = 1.0




Vorwartselimination in python

def elimination(a,b)

n = a.shape [0]

for i in range(n):

for j in range(i+1,n):

faktor = a[j,i]/a[i,i]

b[j] -= b[i]* faktor

for k in range(i,n):

a[j,k] -= a[i,k]* faktor

• 1. Schleife: Uber alle Zeilen/Diagonalelemente• 2. Schleife: Uber alle Zeilen unterhalb des Diagonalelements.

Elemente unterhalb des Diagonalelements mussen Null werden• 3. Schleife: Uber alle Elemente der Zeile aus der 2. Schleife

Vielfaches der Diagonalelementszeile zur aktuellen Zeileaddieren

• faktor so gewahlt, dass Element unterhalb desDiagonalelements zu Null wird




Ruckwartssubstitution in python

def substitution(a,b):

n = a.shape [0]

x = empty(n)

for i in range(n-1, -1, -1):

zaehler = b[i]

for j in range(i+1,n):

zaehler -= a[i,j] * x[j]

x[i] = zaehler / a[i,i]

return x

• 1. Schleife: Zeilen von unten nach oben durchlaufen• 2. Schleife: Bekannte Variablen auf die rechte Seite• Zum Schluss: Teilen durch den Koeffizienten des aktuellen

x-Werts• Losungsvektor x zuruckgeben




Losung des vorigen Beispiels

>>> from lgs import * # enthaelt beide funktionen

# + numpy -import

>>> a = array ([[2,3,1],[4,7,3],[2,4,4]], dtype=float)

>>> b = array ([1,2,2], dtype = float)

>>> elimination(a,b); a; b

array ([[ 2., 3., 1.],

[ 0., 1., 1.],

[ 0., 0., 2.]])

array([ 1., 0., 1.])

>>> substitution(a,b)

array([ 1. , -0.5, 0.5])




Losung der Warmeleitungsgleichung 1D

• Stab mit neun Diskretisierungspunkten• Links beheizt: T = 1, rechter Rand: T = 0

>>> def getwaermelgs1D(n):

>>> a = identity(n)

>>> a = a*-2

>>> for i in range(n-1):

>>> a[i, i+1] = 1

>>> a[i+1, i] = 1

>>> b = zeros(n)

>>> b[0] = -1

>>> return(a,b)

>>>

>>> (a,b) = getwaermelgs1D (9)

>>> loese(a,b)

array ([0.9 , 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1])




Losung der Warmeleitungsgleichung 2D• LGS aus der vorigen Vorlesung• Diskretisierung des Raums mit 3x3 Punkten• Unten beheizt: T = 1, restlicher Rand: T = 0

>>> def getwaermelgs2D ():

>>> a = array([[-4, 1, 0, 1, 0, 0, 0, 0, 0],

>>> [ 1,-4, 1, 0, 1, 0, 0, 0, 0],

...

>>> [ 0, 0, 0, 0, 0, 1, 0, 1,-4]],

>>> dtype = float)

>>> b = array([-1,-1,-1, 0, 0, 0, 0, 0, 0],

>>> dtype = float)

>>> return (a,b)

>>>

>>> (a,b) = getwaermelgs2D ()

>>> loese(a,b)

array ([0.43 , 0.53, 0.43, 0.19, 0.25, 0.19, 0.07, 0.10,

0.07])




Kosten des Gauß-AlgorithmusVorwartselimination• n sei die Anzahl an Unbekannten• Drei verschachtelte Schleifen• Jeweils O(n) Schleifendurchlaufe• Insgesamt Kosten von O(n3)

Ruckwartssubstitution• Zwei verschachtelte Schleifen (je O(n))• Insgesamt Kosten von O(n2)

Optimierung• Matrix ist dunn besetzt, viele der Operationen unnotig• 1D: Matrix tridiagonal

• Beide inneren Schleifen nur uber zwei Elemente• ⇒ Gesamtkosten O(n)

• 2D: Matrix hat zwei weitere Bander mit Abstand√

n• Inneren Schleifen uber je

√n Elemente

• ⇒ Gesamtkosten O(n ·√

n ·√

n) = O(n2)




Iterative VerfahrenMotivation• O(n3) ist viel zu teuer!• Optimiert immer noch O(n2) (2D) bzw. O(n2,33) (3D)• Kosten bei 3D-Simulation:

Mehrgitter Jacobi Gauß opt. GaußN Punkte p. Dim. n = N3 N5 n2,33 = N7 n3 = N9

10 1.000 1e5 1e7 1e920 8.000 3e6 1e9 5e1150 125.000 3e8 8e11 2e15100 1e6 1e10 1e14 1e18

• Idee: Schrittweise Verbesserung einer Naherungslosung xi

• x (k+1) = Φ(x (k))

• Bei Konvergenz gibt es einen Fixpunkt Φ(x) = x .• Ein Schritt des Verfahrens kostet O(n) .• Anzahl Schritte moglichst klein, im Ideallfall O(1)




Typen von IterationsverfahrenRelaxationsverfahren• Ausgangspunkt LGS Ax = b• Ziel: Fehler e(k) = x (k) − x minimieren• Fehler ist leider nicht messbar, da x unbekannt• Residuum r (k) = b − Ax (k) = Ae(k) als Richtung des Fehlers• Neues x unter Verwendung des Residuums• Z.B. bei Jacobi: Residuum lokal immer auf Null setzen

Minimierungsverfahren• Ausgang wieder LGS Ax = b, A symmetrisch, positiv definit• Losung aquivalent zur Minimierung von f (x) := 1

2 xT Ax − bT x + c• f ′(x) := 1

2 AT x + 12 Ax − b = Ax − b = −r(x)

• Minimum von f (Nullstelle von f ′) ist Losung des LGS• Jeder Iterationsschritt lauft einen Schritt (z.B. steilster Abstieg)

auf das Minimum zuT. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird die PDE lokal erfullt.• Neue Werte T (k+1) aus alten Werten T (k) berechnen• 1D Warmeleitung stationar: Zweite Ableitung muss Null sein• ⇒ Ti+1 − 2Ti + Ti−1 = 0 ⇒ Ti = 1

2 (Ti+1 + Ti−1)

• ⇒ T (k+1)i = 1

2 (T (k)i+1 + T (k)

i−1)

T

x

T1

T2

T3

T0

T5

T6

T7

T4




Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird die PDE lokal erfullt.• Neue Werte T (k+1) aus alten Werten T (k) berechnen• 1D Warmeleitung stationar: Zweite Ableitung muss Null sein• ⇒ Ti+1 − 2Ti + Ti−1 = 0 ⇒ Ti = 1

2 (Ti+1 + Ti−1)

• ⇒ T (k+1)i = 1

2 (T (k)i+1 + T (k)

i−1)

T

x

T1

T2

T3

T0

T5

T6

T7

T4




Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird die PDE lokal erfullt• Neue Werte T (k+1) aus alten Werten T (k) berechnen• 1D Warmeleitung stationar: Zweite Ableitung muss Null sein• ⇒ Ti+1 − 2Ti + Ti−1 = 0 ⇒ Ti = 1

2 (Ti+1 + Ti−1)

• ⇒ T (k+1)i = 1

2 (T (k)i+1 + T (k)

i−1)

T

x

T1

T2

T3

T0

T5

T6

T7

T4




Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird die PDE lokal erfullt• Neue Werte T (k+1) aus alten Werte T (k) berechnen• 1D Warmeleitung stationar: Zweite Ableigung muss Null sein• ⇒ Ti+1 − 2Ti + Ti−1 = 0 ⇒ Ti = 1

2 (Ti+1 + Ti−1)

• ⇒ T (k+1)i = 1

2 (T (k)i+1 + T (k)

i−1)

T

x

T1

T2

T3

T0

T5

T6

T7

T4




Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird die PDE lokal erfullt.• Neue Werte T (k+1) aus alten Werten T (k) berechnen• 1D Warmeleitung stationar: Zweite Ableigung muss Null sein• ⇒ Ti+1 − 2Ti + Ti−1 = 0 ⇒ Ti = 1

2 (Ti+1 + Ti−1)

• ⇒ T (k+1)i = 1

2 (T (k)i+1 + T (k)

i−1)

T

x

T1

T2

T3

T0

T5

T6

T7

T4




Jacobi-Verfahren (aus LGS-Sicht)• Matrix A aufteilen

A · x = (M + N) · x = b

• Umformenx = N−1(b −M · x)

• Vorige Gleichung als Iterationsvorschrift

x (k+1) = N−1(b −M · x (k))

= N−1(b − (A− N) · x (k))

= N−1(b − A · x (k) + N · x (k))

= N−1(b − A · x (k)) + x (k)

= x (k) + N−1 · r (k)

• N leicht invertierbar: Matrix der Diagonalelemente von A• x (k+1)

i = x (k)i − 1

2 (−(1 · x (k)i−1 − 2 · x (k)

i + 1 · x (k)i+1)) = 1

2 (x (k)i−1 + x (k)

i+1)

• 2D: x (k+1)i,j = 1

4 (x (k)i−1,j + x (k)

i+1,j + x (k)i,j−1 + x (k)

i,j+1)T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Jacobi-Verfahren ohne Aufstellung der Matrix

• Die Matrix A hat O(n2) Elemente• Die meisten Elemente sind Null⇒ Verschwendung von Speicher• Entweder Matrix effizienter speichern, oder gar nicht aufstellen.• Fur jeden neuen x-Wert sind nur jeweils vier alte (2D) notig• Koeffizient ist dabei jeweils eins• Die Matrix ist also uberhaupt nicht notig.

Randbehandlung

• Die Iterationsvorschrift x (k+1)i,j = ... geht so nur fur innere Punkte.

• Am Rand fehlen Nachbarelemente, dafur muss b berucksichtigtwerden.

• Andere Moglichkeit: Randpunkte in x aufnehmen und nur allebisherigen x-Werte aktualisieren

• Damit zusatzlich zur Matrix A auch noch b unnotig.




Implementierung – Jacobi fur 2D WarmeMaß fur den Fehler• Perfekt ware ein Fehler (damit auch Residuum) gleich Null.• Dazu ist eine sehr hohe Zahl an Iterationen notig.• Eine so hohe Genauigkeit ist meistens nicht erforderlich.• So lange iterieren, bis Residuum unter einer Schranke ist.

• res =√∑n

i=0(resi )2/n

Includes

from numpy import *


from math import sqrt

class Waermesimulation:

...




Implementierung – Jacobi fur 2D WarmeKonstruktor• Vorgabe von Gebietsbreite, -hohe und Zellbreite h• Vorgabe der Temperatur an den vier Randern• Erstellung eines initialen Arrays fur die Temperatur• Null fur innere Zellen, am Rand die gegebenen Werte

def __init__(self , lengthx , lengthy , h, leftT , \

rightT , upperT , lowerT ):

self.zeilen = int(lengthy/h)

self.spalten = int(lengthx/h)

self.T = zeros ((self.zeilen+2,self.spalten +2))

for i in range(1, self.zeilen +1):

self.T[i, 0] = leftT

self.T[i, self.spalten +1] = rightT

for i in range(1, self.spalten +1):

self.T[0, i] = lowerT

self.T[self.zeilen+1, i] = upperT




Iterationsschleife• Iterieren bis Residuum kleiner als maxres

• Nicht direkt auf self.T arbeiten, sondern auf einer Kopie

def jacobi(self , maxres ):

res = maxres +1.

while res > maxres:

tempT = copy(self.T)

res = 0


for j in range(1, self.spalten +1):

res += (tempT[i-1,j] + tempT[i+1,j] \

+ tempT[i,j-1] + tempT[i,j+1] \

- 4 * tempT[i,j])**2

self.T[i,j] = 0.25*( tempT[i-1,j] \

+ tempT[i+1,j] \

+ tempT[i,j-1] \

+ tempT[i,j+1])

res = sqrt(res/(self.zeilen*self.spalten ))




Ergebnis plotten

def plot(self):

print self.T

plt.pcolormesh(self.T)

plt.show()

Beispiel

>>> from waermesim import *

>>> lgs = Waermesimulation (50., 30., 1., 0.6, 0., 1., 0.)

>>> lgs.jacobi (0.01)

>>> lgs.plot()




Gauß-Seidel-Verfahren• Wie Jacobi, es werden aber alle neuen Werte gleich verwendet• Jacobi 1D:⇒ T (k+1)

i = 12 (T (k)

i+1 + T (k)i−1)

• Gauß-Seidel 1D:⇒ T (k+1)i = 1

2 (T (k)i+1 + T (k+1)

i−1 )• Achtung: Gauß-Seidel abhangig von Durchlaufreihenfolge• LGS-Sicht: N = DA + LA (Diagonale + unteres Dreieck)

x (k+1) = x (k) + N−1 · r (k)

T

x

T1

T2

T3

T0

T5

T6

T7

T4





i = 12 (T (k)

i+1 + T (k)i−1)


2 (T (k)i+1 + T (k+1)


x (k+1) = x (k) + N−1 · r (k)

T

x

T1

T2

T3

T0

T5

T6

T7

T4





i = 12 (T (k)

i+1 + T (k)i−1)


2 (T (k)i+1 + T (k+1)


x (k+1) = x (k) + N−1 · r (k)

T

x

T1

T2

T3

T0

T5

T7

T4

T6




Implementierung Gauß-Seidel• Unterschied zu Jacobi: direkt auf self.T arbeiten• Vorteil: schneller und speichereffizienter

def gaussseidel(self , maxres ):

res = abs(maxres )+1.

while abs(res) > abs(maxres ):

#tempT = copy(self.T)

res = 0



res += (self.T[i-1,j] + self.T[i+1,j] \

+ self.T[i,j-1] + self.T[i,j+1] \

- 4 * self.T[i,j])**2

self.T[i,j] = 0.25*( self.T[i-1,j] \

+ self.T[i+1,j] \

+ self.T[i,j-1] \

+ self.T[i,j+1])

res = sqrt(res/(self.zeilen*self.spalten ))




Werbung: Chair of Informatics V—SCCS

Efficient Numerical Algorithms—Parallel & HPC

• High-dimensional numerics (sparse grids)• Fast iterative solvers (multi-level methods,

preconditioners)• Adaptive, octree-based grid generation• Space-filling curves• Numerical linear algebra• Numerical algorithms for image processing• HW-aware numerical programming

Fields of application in simulation

• CFD (incl. fluid-structure interaction)• Plasma physics• Molecular dynamics• (Quantum chemistry)Further info→ www5.in.tum.deFeel free to come around and ask for thesis topics ;-)




Teil XV

Mehrgitterverfahren zur Losunglinearer Gleichungssysteme




Wiederholung: Relaxationsverfahren

• Direkte Losung des LGS ist viel zu teuer;• Daher wird versucht, den Fehler lokal zu beseitigen.• Echter Fehler nicht messbar, daher Zugang uber Residuum;• Innere Schleife lauft uber alle Gitterpunkte;• Außere Schleife wiederholt Relaxation bis das Residuum klein

genug ist;• Anzahl der außeren Schleifendurchlaufe von der Anzahl an

Gitterpunkten abhangig;• Auf den folgenden Folien: Nebeneinanderstellung von PDE-Sicht

und LGS-Sicht.




u uuu uuuu u Jacobi-Verfahren:

xnew = xold + D−1 ·(

b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uu uu Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uu u uu Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uu uu

Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uuu uu

Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uuu u uu

Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uuu uu Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uuu u

Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uuu u u

Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




u uuu uuuu uuu u

Jacobi-Verfahren:


b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9

− h2

− 14 0 0 0 0 0 0 0 0

0 − 14 0 0 0 0 0 0 0

0 0 − 14 0 0 0 0 0 0

0 0 0 − 14 0 0 0 0 0

0 0 0 0 − 14 0 0 0 0

0 0 0 0 0 − 14 0 0 0

0 0 0 0 0 0 − 14 0 0

0 0 0 0 0 0 0 − 14 0

0 0 0 0 0 0 0 0 − 14

·

f1f2f3f4f5f6f7f8f9

−1

h2

−4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

·

uold1

uold2

uold3

uold4

uold5

uold6

uold7

uold8

uold9




Jacobi-Verfahren: Interaktive Rechnung

u u uT0 T0 T0

u uuu uuuu uuu u

1

4

7

2

5

8

3

6

9Jacobi-Verfahren:

xnew = xold +D−1·(

b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

48840

−4−8−8−4

+1

4

161616000000

+1

4

1 2 3 4 5 6 7 8 9

− 4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

1 2 3 4 5 6 7 8 9

·

4 18 28 34 40 5

−4 6−8 7−8 8−4 9

unew =

1 2 3 4 5 6 7 8 9

7 7 5 − 1 0 1 − 1 − 3 − 4

T




Jacobi-Verfahren: Interaktive Rechnung – korrekt;-)

u u uT0 T0 T0

u uuu uuuu uuu u

1

4

7

2

5

8

3

6

9Jacobi-Verfahren:

xnew = xold +D−1·(

b− A · xold)

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

48840

−4−8−8−4

+1

4

161616000000

+1

4

1 2 3 4 5 6 7 8 9

− 4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

1 2 3 4 5 6 7 8 9

·

4 18 28 34 40 5

−4 6−8 7−8 8−4 9

unew =

1 2 3 4 5 6 7 8 9

7 7 5 −1 0 1 −1 −3 −3

T




Jacobi-Verfahren: Interaktive Rechnung II

u u uT0 T0 T0

u uuu uuuu uuu u

1

4

7

2

5

8

3

6

9

u u u16 16 16u uuu u

uuu uuu u

7

-1

-1

7

0

-3

5

1

-3

Jacobi-Verfahren:

xnew = xold +D−1·

b− A · xold︸︷︷︸=resi=−Sterni

xnewi = xold

i −1aii

[Stern]i

unew1

unew2

unew3

unew4

unew5

unew6

unew7

unew8

unew9

=

775

−101

−1−3−3

+1

4

161616000000

+1

4

1 2 3 4 5 6 7 8 9

− 4 1 0 1 0 0 0 0 01 −4 1 0 1 0 0 0 00 1 −4 0 0 1 0 0 01 0 0 −4 1 0 1 0 00 1 0 1 −4 1 0 1 00 0 1 0 1 −4 0 0 10 0 0 1 0 0 −4 1 00 0 0 0 1 0 1 −4 10 0 0 0 0 1 0 1 −4

1 2 3 4 5 6 7 8 9

·

7 17 25 3

−1 40 51 6

−1 7−3 8−3 9

unew =

1 2 3 4 5 6 7 8 9

5.5 7 6 1.5 1 0.5 − 1 − 1 − 0.5

T

;-)T. Neckel: Einfuhrung in die wissenschaftliche Programmierung



Konvergenz Jacobi: Fourier-Analyse

Szenario• Stationare Warmeleitung bei quadratischer Flache• Randbedingung uberall Null• Losung: Null auf der gesamten Flache• Initialer Losungsvektor: x0 = e0

• Untersuchung: Veranderung des Fehlers in einer Iteration• Initialer Fehler als Kombination von Sinus-Kurven

unterschiedlicher Frequenz• ab jetzt zur Vereinfachung: 1D: x =

∑m am · (sin(mπhi))i=1..N




Was sind hohe/niedrige Frequenzen?

• Niedrige Frequenz: eine Schwingung (sin(πx)) im betrachtetenBereich

• Ist sin(1000πx) eine hohe Frequenz?

Nyquist-Shannon-Abtasttheorem• Eine Frequenz ist hoch/niedrig bezuglich der Anzahl an

Abtastpunkten.• Ein Signal mit Maximalfrequenz fmax muss mit einer Frequenz

großer als 2 · fmax abgetastet werden.• Fur gegebene Abtastfrequenz f ist die hochste darstellbare

Frequenz also 1/2f .• Die Abtastfrequenz entspricht gerade der Anzahl an

Diskretisierungspunkten!




Anderung des Fehlers in einer Iteration• x (k+1) = x (k) − D−1Ax (k) = (I − D−1A)x (k) = Mx (k)

• e(k+1) = Me(k)

• (M(sin(mπhi)))i = cos(mπh)(sin(mπhi))i

Eigenwerte von M• λm = cos(mπh) sind Eigenwerte zu den Eigenvektoren

(sin(mπhi))i• Frequenzen zu betragsmaßig niedrigen Eigenwerten (nahe Null)

verschwinden schnell.• Frequenzen zu betragsmaßig hohen Eigenwerten (nahe Eins)

verschwinden langsam.• Fur x gegen Null nimmt cos(x) den maximalen Wert an• Taylor-Entwicklung um 0 + πh:

λ1 = cos(0)− πhsin(0)− cos(0)π2h2/2 + ...

≈ 1− 0− π2h2/2= 1− c · h2




Konvergenz der Iteration

Was passiert bei Verdopplung der Gitterpunkte?• Eine Iteration mit h druckt den Fehler um 1− c · h2

• h→ h/2• Eine Iteration druckt Fehler um: 1− c · (h/2)2 = 1− (1/4)c · h2

• Zwei Iterationen: (1− c · (h/2)2)2 ≈ 1− (1/2)c · h2

• Vier Iterationen: (1− (1/2)c · h2)2 ≈ 1− c · h2

• Bei Halbierung von h sind vier Mal so viele Iterationen notig!• ⇒ Anzahl an Iterationen wachst mit 1/h2

• Wie konnen wir es besser machen?




Demo “Militarkapelle”




Reduktion des Fehlers bei Jacobi/Gauss-Seidel• Jacobi-Verfahren erfullt die PDE lokal;• D.h. lokal ist die zweite Ableitung Null bzw. sehr klein.• D.h. fur glatte Funktionen ist das Residuum gering.• Fur hohe Frequenzen ist das Residuum groß, d.h. solche

Frequenzen werden schnell beseitigt, niedrige nicht.• Abbildung: 1D Warmeleitung mit Rand=0 und zufalligen

Startwerden fur die Temperatur

1

-0.75

0

0.75

1

-0.75

0

0.75

Jacobi, 10 SchritteJacobi, 100 SchritteMehrgitter, 1 Schritt




Mehrgitter: Motivation

• Optimales LGS-Verfahren:• Rechenaufwand wachst linear mit der Anzahl an

Diskretisierungspunkten• Besser geht es nicht (jeder Punkt muss mindestens einmal

betrachtet/bearbeitet werden)• Jacobi/Gauss-Seidel:

• Kosten einer Iteration bei sind O(n)• Die Anzahl an Iterationen wachst mit der Anzahl an

Diskretisierungspunkten.• Niedrige Fehlerfrequenzen verschwinden nur sehr langsam.• Ursache: Informationen mussen durch das gesamte Gebiet

(damit durch alle Punkte) propagiert werden• Mehrgitter-Verfahren propagiert Informationen sehr viel schneller.




Grundidee des Mehrgitter-Verfahrens

• Jacobi/Gauss-Seidel sind schlecht fur niedrige Frequenzen• Shannon: Eine Frequenz ist niedrig bezuglich einer bestimmten

Zahl an Diskretisierungspunkten• Bei Halbierung der Diskretisierungspunkte wird die Frequenz

bezuglich der Diskretisierung verdoppelt• D.h. fur jede Frequenz gibt es eine Diskretisierung, fur die diese

Frequenz mittel/hoch ist• Durch Einsatz verschiedener Gitterauflosungen gibt es fur jede

Komponente des Fehlers ein Gitter, bezuglich dessen dieFrequenz der Fehlerkomponente niedrig ist.




Grundidee des Mehrgitter-Verfahrens

Tafelbild




Mehrgitter-Algorithmus (generell, Pseudocode)

def mehrgitter(level , b, x):

# Vorglaetten (hohe Fehlerfrequenzen beseitigen)

self.glaetter(level , b, x)

if(level < maxlevel ):

# Berechnung des Residuums

r = residuum(level , b, x)

# Uebertragung des Residuums auf das grobe Gitter

b_c = restriktion(level , r)

# Loesung von -Ae = r auf grobem Gitter

e_c = mehrgitter(level+1, b_c , zero_vector(level))

# Prolongation

x_delta = prolongation(level , e_c)

# Korrektur

x = x + x_delta

# Nachglaetten / oberster level

self.glaetter(2, level)




Restriktion

• Fur neue Naherung xk+1 auf grobem Gitter sind notig: xk und r k

bezuglich des groben Gitters• Bisherige Naherungslosung muss auf das nachstgrobere Gitter

transportiert werden• Einfachste Methode: Injektion. Nur jeder zweite Punkt wird

verwendet, die dazwischen verworfen• Teilweise bessere Ergebnisse, wenn zwischen drei

benachbarten Punkten gemittelt wird




Prolongation

• Der auf dem groben Gitter berechnete Fehler muss auf das feineubertragen werden

• Feingitterpunkte, die zugleich Grobgitterpunkte sind, erhaltenden Wert direkt vom Grobgitterpunkt

• Feingitterpunkte zwischen zwei Grobgitterpunkte erhalten denWert durch Interpolation




Mehrgitter-Algorithmus (vereinfacht)• Im bisherigen Code (letzte Vorlesung) zur Warmeleitung kein b

vorgesehen• Daher nicht Losung von −Ae = r auf grobem Gitter, sondern

direkte Losung

def mehrgitter_rekursiv(self , level):

if(2** level < min(self.zeilen , self.spalten )):

# Vorglaetten


# Restriktion: nicht notig (Injektion)

# -

# Loesung auf grobem Gitter

self.mehrgitter_rekursiv(level +1)

# Prolongation

self.prolongation(level)

# Nachglaetten / oberster Level





Glaetter und Prolongation

def glaetter(self , n, l):

offset = 2**l

for iters in range(n):

for i in range (2**l, self.zeilen+1, 2**l):

for j in range (2**l, self.spalten+1, 2**l):

self.T[i,j] = 0.25*( self.T[i-offset ,j] \

+ self.T[i+offset ,j] \

+ self.T[i,j-offset] \

+ self.T[i,j+offset ])

def prolongation(self , l):

offset = 2**l

for i in range (2**l, self.zeilen+1, 2**(l+1)):

for j in range (2**l, self.spalten+1, 2**(l+1)):

self.T[i,j] = 0.25*( self.T[i-offset ,j] \

+ self.T[i+offset ,j] \

+ self.T[i,j-offset] \

+ self.T[i,j+offset ])




Rekursionsaufruf

def mehrgitter(self , maxres ):

res = self.getResiduum ()

while res > maxres:

self.mehrgitter_rekursiv (0)

res = self.getResiduum ()

Berechnung des Residuum

def getResiduum(self):

res = 0



res += (self.T[i-1,j] + self.T[i+1,j] \

+ self.T[i,j-1] + self.T[i,j+1] \

- 4 * self.T[i,j])**2

return sqrt(res/(self.zeilen*self.spalten ))




Beispiel

>>> from waermesim import *

>>> lgs = Waermesimulation (4 , 5 ,1. ,0.6 ,0. ,1. ,0.)

>>> lgs.mehrgitter (0.001)

>>> lgs.plot()




Kosten des Mehrgitter-Verfahrens (2D)

• Pro Level zwei bis vier Jacobi/Gauss-Seidel Iterationen: O(nl )

• Residuumsberechnung: O(nl )

• Restriktion: O(nl )

• Prolongation: O(nl )

Rekursion• Kosten auf feinstem Level: cn• Kosten auf zweitfeinstem Level: 1/4cn• ...• ⇒ Kosten cn

∑∞i=0(1/4)i = 4/3cn

Anzahl Iterationen• O(1)

• ⇒ Gesamtkosten: O(n)




KlausurstoffTeil I: Erste Schritte jaTeil II: Datentypen jaTeil III: Kontrollstrukturen jaTeil IV: Funktionen und Module jaTeil V: IO und Datentypen ja (ausser Modul os)Teil VI: OOP ja (außer UML)Teil VII: Regulare Ausdrucke Theorie nicht, Verwendung schon

(Maskierungszeichen wie \d, ... werdenangegeben)

Teil VIII: Exceptions ja (außer Build-In Exceptions)Teil IX: Grafik Nein (außer Rekursion)Teil X: Partikelsysteme jaTeil XI: Baume, ... jaTeil XII: Wiss. Rechnen Nur numpy, keine PaketnamenTeil XIII: Warmeleitung neinTeil XIV: Losung von LGS jaTeil XV: Mehrgitterverfahren nein



IN8008 Organisatorisches und erste Wintersemester 2013 ... · PDF fileIN8008, Wintersemester...

Documents

Transcript of IN8008 Organisatorisches und erste Wintersemester 2013 ... · PDF fileIN8008, Wintersemester...