Einfuhrung¨ in die wissenschaftliche Programmierung · PDF fileScientific Computing in...
Transcript of Einfuhrung¨ in die wissenschaftliche Programmierung · PDF fileScientific Computing in...
Scientific Computing in Computer Science, Technische Universitat Munchen
Einfuhrung in die wissenschaftlicheProgrammierung
IN8008
Tobias Neckel
Wintersemester 2011/2012
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 1
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil I
Organisatorisches und erste Schritte
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 2
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 3
Scientific Computing in Computer Science, Technische Universitat Munchen
Organisatorisches
Raum und Zeit• Modul IN8008• Vorlesung Mo 10:10 - 11:30, MI HS1• Zentralubung Mi 08:30 - 10:00, PH HS1
Informationen• Webseite: www5.in.tum.de→ Teaching→ Einfuhrung in die
wissenschaftliche Programmierung• FRAGEN!
Personen• Dr. Andreas Paul ([email protected]): Ubungsleitung• Benjamin Peherstorfer ([email protected]): Ubungsleitung• Dr. Tobias Neckel ([email protected]): Vorlesung
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 4
Scientific Computing in Computer Science, Technische Universitat Munchen
Organisatorisches IIBonus• 2 Programmieraufgaben-Blatter (zusatzlich zu
ZU-Hausaufgaben)• Erfolgreiche Bearbeitung + Prasentation• 3er-Gruppen
Notige Software• Python 2.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
ToDo• TUMonline: Anmelden fur Veranstaltung• TUMonline: Anmelden fur Klausur
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 5
Scientific Computing in Computer Science, Technische Universitat Munchen
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, ...)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 6
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 7
Scientific Computing in Computer Science, Technische Universitat Munchen
Warum Python? Warum Python in der Physik?
• Viele Physiker verwenden es (MPI etc.).• Pre- und Postprocessing sehr elegant!• Scripting!
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 8
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 9
Scientific Computing in Computer Science, Technische Universitat Munchen
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
Interactive Hilfe• help() fur interactive Hilfe• help(’modules’) Liste der verfugbaren Module• help(object|’command’) Hilfe zu Objekt oder Befehl
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 10
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 11
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 12
Scientific Computing in Computer Science, Technische Universitat Munchen
Python ausfuhren: python. . .• Python interpreter starten (python oder ipython)• Python-prompt: >>> Befehl eintippen und mit <Enter>
bestatigen• ipython-prompt: In [1]:
$ python
Python 2.5.2 (r252 :60911 , Jan 20 2010, 23:16:55)
[GCC 4.3.2] 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.5! This is the online help
utility. [...]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 13
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 14
Scientific Computing in Computer Science, Technische Universitat Munchen
Python ausfuhren: Datei ubergeben
• Datei square_me.py:
print 6*2
• Ausfuhren mit
$ python square_me.py
12
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 15
Scientific Computing in Computer Science, Technische Universitat Munchen
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
12
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 16
Scientific Computing in Computer Science, Technische Universitat Munchen
“Hello World!”Python – lesbarer Code• In Python
>>> print "Hello World!"
Hello World!
• In Java
public class Hello
public static void main(String argv [])
System.out.println("Hello World!");
$ javac Hello.java
$ java Hello
Hello World!
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 17
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 18
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 19
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 20
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil II
Datentypen
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 21
Scientific Computing in Computer Science, Technische Universitat Munchen
Numerische Typen - ganze Zahlen• 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’>
• Achtung beim Dividieren: Nachkommastellen werden abgehackt
>>> 5/3
1
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 22
Scientific Computing in Computer Science, Technische Universitat Munchen
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”)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 23
Scientific Computing in Computer Science, Technische Universitat Munchen
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• Alles in Python ist ein Objekt• Zugriff auf Objektmethoden: <Objektname>.<Methodenname>
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 24
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 25
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 26
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 27
Scientific Computing in Computer Science, Technische Universitat Munchen
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")
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 28
Scientific Computing in Computer Science, Technische Universitat Munchen
Ausgabe verschiedener Datentypen
>>> 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 29
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 30
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 31
Scientific Computing in Computer Science, Technische Universitat Munchen
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]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 32
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 33
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil III
Kontrollstrukturen
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 34
Scientific Computing in Computer Science, Technische Universitat Munchen
Struktogramme
• Graphische Darstellung von Programmen• Anweisungen
Anweisung 1Anweisung 2
...Anweisung n
• BedingungenBedingung
wahr falschAnweisungs-
block 1Anweisungs-
block 2
• SchleifenSchleifenbedinung
Anweisungsblock
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 35
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 36
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 37
Scientific Computing in Computer Science, Technische Universitat Munchen
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’
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 38
Scientific Computing in Computer Science, Technische Universitat Munchen
Nochmal Typen und Objekte
• 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’>
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 39
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 40
Scientific Computing in Computer Science, Technische Universitat Munchen
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; exp = 0
while a > 1:
a /= 2 # a = a / 2
exp+=1 # exp = exp + 1
print a, "=", 2, "^", exp
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 41
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 42
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Numerische Losung von ODEs
• Logistische Differentialgleichung: p(t) = ap(t)− bp(t)2
• Anfangswert p(0) = p0
• Losung p(t) = a·p0b·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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 43
Scientific Computing in Computer Science, Technische Universitat Munchen
Mittels 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?
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 44
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 45
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 46
Scientific Computing in Computer Science, Technische Universitat Munchen
Das Programm fur das Euler-Verfahren sieht ganz ahnlich aus wiedas fur die Wertetabelle von vorhin:
a = 1. # Parameter der DGL
b = 0.01
p0 = 10. # Anfangswert
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 47
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 48
Scientific Computing in Computer Science, Technische Universitat Munchen
Erweitertes Programm
>from math import exp
a = 1. # Parameter der DGL
b = 0.01
p0 = 10. # Anfangswert
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 49
Scientific Computing in Computer Science, Technische Universitat Munchen
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 46
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 50
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 51
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 52
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 53
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 54
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil IV
Funktionen und Module
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 55
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ...)• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 56
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ...)• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 56
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 57
Scientific Computing in Computer Science, Technische Universitat Munchen
Funktionen sind Objekte
>>> def faculty(n):
>>> fac = 1
>>> for i in range(2, n+1):
>>> fac = fac*i
>>> return fac
>>>
>>> faculty (3)
6
>>> f = faculty
>>> f(4)
24
>>> type(f)
<type ’function ’>
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 58
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 59
Scientific Computing in Computer Science, Technische Universitat Munchen
Namensparameter• Bisher: Folge von Aktualparametern gemaß Reihenfolge der
Formalparameter (Positionsparameter)• Durch Angabe von Folge Formalparameter=Aktualparameter bei
Aufruf beliebige Reihenfolge moglich (Namensparameter)
>>> x = potenziere(basis=2, exponent =10)
>>> x = potenziere(exponent =10, basis =2)
• Mischung von Positions- und Namensparameter moglich• Erst Folge von Positions- dann die verbleibenden
Namensparameter
>>> x = potenziere (2, exponent =10)
>>> x = potenziere(basis=2, 10) # Fehler
>>> x = potenziere (2, basis =2) # Fehler
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 60
Scientific Computing in Computer Science, Technische Universitat Munchen
Standardwerte fur Parameter• Namensparameter relativ nutzlos, wenn alle angegeben werden
mussen• 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 kann es hilfreich sein
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 61
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 62
Scientific Computing in Computer Science, Technische Universitat Munchen
GultigkeitsbereichFunktionen• 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 außerhalb nicht sichtbar.• Grundregel: GLOBALE VARIABLEN VERMEIDEN!!!• Variablen, die in einer Funktion definiert oder links in einer
Zuweisung verwendet werden, sind lokal.• Globale Variablen konnen gelesen werden.• Um globale Variablen in einer Funktion zu schreiben, wird das
Schusselwort global verwendet.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 63
Scientific Computing in Computer Science, Technische Universitat Munchen
GultigkeitsbereichFunktionen• 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 außerhalb nicht sichtbar.• Grundregel: GLOBALE VARIABLEN VERMEIDEN!!!• Variablen, die in einer Funktion definiert oder links in einer
Zuweisung verwendet werden, sind lokal.• Globale Variablen konnen gelesen werden.• Um globale Variablen in einer Funktion zu schreiben, wird das
Schusselwort global verwendet.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 63
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiele
>>> def zaehle(n):
>>> global anzahl
>>> anzahl += 1
>>>
>>> anzahl = 0
>>> zaehle (); zaehle (); zaehle ()
>>> anzahl
3
>>> def no_effect ():
>>> xyz = global **2 # Nochmal: ganz schlechter Stil
>>>
>>> global = 9
>>> no_effect ()
>>> xyz
NameError: name ’xyz’ is not defined
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 64
Scientific Computing in Computer Science, Technische Universitat Munchen
Docstrings• Ein Grund fur Funktionen: Funktionalitat kapseln• Jemand, der die Funktion nicht geschrieben hat, soll sie
verwenden konnen.• Dazu ist Dokumentation notig.• Ein String (ein- oder mehrzeilig) nach dem Header ist dafur
vorgesehen.• Python ignoriert den String (Ausdruck ohne Zuweisung).• 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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 65
Scientific Computing in Computer Science, Technische Universitat Munchen
Call by Reference vs. Call by value
• Call by Reference: ”Adresse“ (id) wird ubergeben• Call by Value: Wert der Variablen wird ubergeben• In Python immer Call by Reference• Allerdings: Zuweisung andert nur Referenz innerhalb der
Funktion
>>> def add_one(x):
>>> x = x + 1
>>>
>>> x = 3
>>> add_one(x)
>>> x
3
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 66
Scientific Computing in Computer Science, Technische Universitat Munchen
Call by Reference vs. Call by value (2)
• Nach der Zuweisung zeigt die lokale Variable auf ein neuesObjekt.
• Mit Methoden kann das bisherige Objekt 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]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 67
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 68
Scientific Computing in Computer Science, Technische Universitat Munchen
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]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 69
Scientific Computing in Computer Science, Technische Universitat Munchen
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]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 69
Scientific Computing in Computer Science, Technische Universitat Munchen
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)]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 70
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 71
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Numerische Integration (Quadratur)Berechnung des Integrals
F1(f ,a,b) :=
∫ b
af (x) dx
Rechtecksregel• Bestimmung des Funktionswertes in der Intervallmitte• Multiplikation mit der Intervallbreite• F1 ≈ R := (b − a) · f ( a+b
2 )
• Polynom 0. Ordnung
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 72
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Numerische Integration (Quadratur)Berechnung des Integrals
F1(f ,a,b) :=
∫ b
af (x) dx
Rechtecksregel• Bestimmung des Funktionswertes in der Intervallmitte• Multiplikation mit der Intervallbreite• F1 ≈ R := (b − a) · f ( a+b
2 )
• Polynom 0. OrdnungTrapezregel• 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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 72
Scientific Computing in Computer Science, Technische Universitat Munchen
Simpsonregel• Anderer Name: Keplersche Fassregel• Bestimmung des Funktionswertes an beiden Randern und in der
Mitte• Polynom 2. Ordnung durch diese drei Punkte• Flache unter dem Polynom wird berechnet
• F1 ≈ S := (b − a) · f (a)+4·f ( a+b2 )+f (b)
6
Quadraturfehler
|T − F1| ≤112· sup
x∈[a,b]|f ′′(x)| · (b − a)3
|S − F1| ≤1
2880· sup
x∈[a,b]|f (4)(x)| · (b − a)5
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 73
Scientific Computing in Computer Science, Technische Universitat Munchen
Simpsonregel• Anderer Name: Keplersche Fassregel• Bestimmung des Funktionswertes an beiden Randern und in der
Mitte• Polynom 2. Ordnung durch diese drei Punkte• Flache unter dem Polynom wird berechnet
• F1 ≈ S := (b − a) · f (a)+4·f ( a+b2 )+f (b)
6
Quadraturfehler
|T − F1| ≤112· sup
x∈[a,b]|f ′′(x)| · (b − a)3
|S − F1| ≤1
2880· sup
x∈[a,b]|f (4)(x)| · (b − a)5
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 73
Scientific Computing in Computer Science, Technische Universitat Munchen
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
]
Simpsonsumme
SS :=h3
[f (a) + 4f (a + h) + 2f (a + 2h) + 4f (a + 3h) + . . .
+ 4f (b − h) + f (b)]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 74
Scientific Computing in Computer Science, Technische Universitat Munchen
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
]
Simpsonsumme
SS :=h3
[f (a) + 4f (a + h) + 2f (a + 2h) + 4f (a + 3h) + . . .
+ 4f (b − h) + f (b)]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 74
Scientific Computing in Computer Science, Technische Universitat Munchen
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
]
Simpsonsumme
SS :=h3
[f (a) + 4f (a + h) + 2f (a + 2h) + 4f (a + 3h) + . . .
+ 4f (b − h) + f (b)]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 74
Scientific Computing in Computer Science, Technische Universitat Munchen
Fehler der Trapezsumme• 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
Fehler der Simpsonsumme• Analog zur Trapessumme ergibt sich
|SS − F1| ≤1
2880· sup
x∈[a,b]|f (4)(x)| · (b − a) · h4
• Verdopplung des Rechenaufwands (Halbierung von h) reduziertden maximalen Fehler um den Faktor 16
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 75
Scientific Computing in Computer Science, Technische Universitat Munchen
Fehler der Trapezsumme• 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
Fehler der Simpsonsumme• Analog zur Trapessumme ergibt sich
|SS − F1| ≤1
2880· sup
x∈[a,b]|f (4)(x)| · (b − a) · h4
• Verdopplung des Rechenaufwands (Halbierung von h) reduziertden maximalen Fehler um den Faktor 16
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 75
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 76
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 76
Scientific Computing in Computer Science, Technische Universitat Munchen
Berechnung mit unterschiedlicher Genauigkeit• 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 77
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 78
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 79
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 80
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 80
Scientific Computing in Computer Science, Technische Universitat Munchen
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"
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 81
Scientific Computing in Computer Science, Technische Universitat Munchen
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"
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 81
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 82
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 82
Scientific Computing in Computer Science, Technische Universitat Munchen
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"]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 83
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil V
IO und weitere Datentypen
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 84
Scientific Computing in Computer Science, Technische Universitat Munchen
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, . . . )
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 85
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 86
Scientific Computing in Computer Science, Technische Universitat Munchen
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 = f.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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 87
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 88
Scientific Computing in Computer Science, Technische Universitat Munchen
Kommandozeilenparameter
sys.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,
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 89
Scientific Computing in Computer Science, Technische Universitat Munchen
Modul optparse• Klasse OptionParser im module optparse vereinfacht das Parsen
von Optionen
#!/usr/bin/python
import optparse
p = optparse.OptionParser ()
# Optionen spezifizieren
p.add_option("-o", action="store", dest="opt")
p.add_option("-v","--verbose",action="store_true",
dest="verbose",
help="Erzeugt ausfuhrliche Ausgabe")
# Optionen parsen
(options , args) = p.parse_args ()
if options.verbose:
print "Programm startet ..."
print options , args
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 90
Scientific Computing in Computer Science, Technische Universitat Munchen
Modul optparse (2)• Das Programm unterstutzt jetzt
• Standardhilfe uber Parameter -h• Zwei vordefinierte Parameter
• Ausprobieren (cmdline.py):
./ cmdline.py
./ cmdline.py -h
./ cmdline.py -o "Blubb" --verbose
./ cmdline.py -v file1.txt file2.txt
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 91
Scientific Computing in Computer Science, Technische Universitat Munchen
Ausgewahlte Parameter von add_option
• Jede Option ("--xyz") kann eine Kurzform ("-x") haben• action: ’store’ (default), ’store_true’, ’store_false’, ’append’,
’count’, ...• dest – Name der Variable zum Speichern der Option• default – Standardwert• help – Hilfetext• type – Typ der Variablen: ’string’ (default), ’int’, ’long’,
’float’, ’complex’, ’choice’• choices = [’first’, ’second’, ’third’] – fur type=’choice’
Hilfetext mit set_usage• Referenz auf den Programmnamen mit %prog
p.set_usage(’’’%prog [optionen] Datei(en)
Macht rein gar nichts mit den Dateien ’’’)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 92
Scientific Computing in Computer Science, Technische Universitat Munchen
nochmal Sequenzen: Dictionaries und Mengen
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 bestimmten Titel
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 93
Scientific Computing in Computer Science, Technische Universitat Munchen
Dictionaries• 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"
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 94
Scientific Computing in Computer Science, Technische Universitat Munchen
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]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 95
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ’
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 96
Scientific Computing in Computer Science, Technische Universitat Munchen
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: -, ˆ, |, &
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 97
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 98
Scientific Computing in Computer Science, Technische Universitat Munchen
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 , ...)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 99
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 100
Scientific Computing in Computer Science, Technische Universitat Munchen
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
...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 101
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 102
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 103
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 104
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 105
Scientific Computing in Computer Science, Technische Universitat Munchen
• 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 ())
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 106
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil VI
Objektorientierte Programmierung
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 107
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 108
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Liste umdrehenAls Funktion
def reverse(liste):
kopie = list(liste)
for i in range(len(liste )):
liste[i] = kopie[len(liste)-i-1]
l = [1, 2, 3]
reverse(l)
Als Methode
Class list()
def reverse(self):
...
l = [1, 2, 3]
l.reverse ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 109
Scientific Computing in Computer Science, Technische Universitat Munchen
Wie programmiert man objektorientiert?
• Andert eure Denkweise• Andert euren Blickpunkt• Es ist keine rein technische Frage
Besonderheiten in python?• Ihr verwendet schon die ganze Zeit objekorientierte Konzepte• Ihr merkt 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’
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 110
Scientific Computing in Computer Science, Technische Universitat Munchen
Theorie: Klassen und Objekte
Klassen:”Idee“ aller Objekte einer Art;Beispiel: Auto
Instanzen / Objekte:Konkretes Ding;Beispiel: ”Das rote Auto dort“
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 111
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 112
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 113
Scientific Computing in Computer Science, Technische Universitat Munchen
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()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 114
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 115
Scientific Computing in Computer Science, Technische Universitat Munchen
Erstes Klassenbeispiel
#!/usr/bin/python
class Student(object ):
anzahl = 0
def __init__(self , name): # Konstruktor
self.name = name # Objektvariable
Student.anzahl += 1 # Klassenvariable
def getName(self): # Methode
return self.name
@staticmethod
def getAnzahl ():
return Student.anzahl
person1 = Student("Alex")
person2 = Student("Michael")
print person1.getName () # "Alex"
print Student.getAnzahl () # 2
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 116
Scientific Computing in Computer Science, Technische Universitat Munchen
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])
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 117
Scientific Computing in Computer Science, Technische Universitat Munchen
Objektorientierter Ansatz: Ofen
#!/usr/bin/python
class Ofen(object ):
def __init__(self , temperatur , modus):
self.temperatur=temperatur
self.modus=modus
def backe ():
Anweisungsblock
meinOfen = Ofen (180 ,2)
meinOfen.backe ()
Und wieder die ganze Backerei
ofenliste = []
ofenliste.append(Ofen (180 ,2))
...
for ofen in ofenliste:
ofen.backe ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 118
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 119
Scientific Computing in Computer Science, Technische Universitat Munchen
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• Syntax: beginnend mit __, aber nicht endend mit __
(__privateVar, __privateMet())• Achtung: indirekt zugreifbar (._Klasse__privateVar)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 120
Scientific Computing in Computer Science, Technische Universitat Munchen
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 +1
def __mul__(self ,b):
return self.value*b.value *2
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 121
Scientific Computing in Computer Science, Technische Universitat Munchen
Vererbung in der Realitat
Objekt
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 122
Scientific Computing in Computer Science, Technische Universitat Munchen
Vererbung in der Realitat
Objekt
Maschine
Tier
Möbel
Pflanze
Comput.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 122
Scientific Computing in Computer Science, Technische Universitat Munchen
Vererbung in der Realitat
Objekt
Maschine
Auto
Tier
Säuget.
Möbel
StuhlTisch
Pflanze
Baum
Blume
Comput.
Laptop
Desktop
Vectorc.
Cluster
Reptilien
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 122
Scientific Computing in Computer Science, Technische Universitat Munchen
Vererbung in der Realitat
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 122
Scientific Computing in Computer Science, Technische Universitat Munchen
Vererbung in der Realitat
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 122
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 123
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 124
Scientific Computing in Computer Science, Technische Universitat Munchen
Strukturdiagramme• Klassendiagramm• Objektdiagramm• Kompositionsstrukturdiagramm• Komponentendiagramm• Verteilungsdiagramm• Paketdiagramm• Profildiagramm
Verhaltensdiagramme• Aktivitatsdiagramm• Anwendungsfalldiagramm (Use-Cases)• Interaktionsubersichtsdiagramm• Kommunikationsdiagramm• Sequenzdiagramm• Zeitverlaufsdiagramm• Zustandsdiagramm
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 125
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 126
Scientific Computing in Computer Science, Technische Universitat Munchen
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
:)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 127
Scientific Computing in Computer Science, Technische Universitat Munchen
Vererbung• Darstellung durch Pfeil mit ”geschlossener“Spitze•
”ist ein“- Beziehung• Pfeil von spezialisierter Klasse zur Vaterklasse
Fahrzeug
Auto Motorrad
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 128
Scientific Computing in Computer Science, Technische Universitat Munchen
Assoziation• Semantischer Zusammenhang zwischen Klassen• Beschreibung der Assoziation neben dem Pfeil
Vorlesung
Dozent Student
liest hört
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 129
Scientific Computing in Computer Science, Technische Universitat Munchen
Aggregation und Komposition• Aggregation (leere Raute): ”hat ein“• Komposition (volle Raute): ”enthalt ein“• Angabe der Anzahl moglich
Fahrzeug
Besitzer Türen1 5
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 130
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil VII
Regulare Ausdrucke
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 131
Scientific Computing in Computer Science, Technische Universitat Munchen
Formale SprachenSyntax• Jede Sprache (naturliche und formale) hat Regeln• Fur einen Text lasst sich feststellen, ob er zu einer Sprache
gehort:• Deutsch: Der Student sitzt 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 dem Studenten.• Python: print "sin(x) = %f"% cos(x)
• Eine Semantik-Korrektur fur Programmiersprachen gibt es nochnicht
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 132
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ><Objt >
AB --> AaB
• S ∈ (V\Σ) ist das Startsymbol
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 133
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 134
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 135
Scientific Computing in Computer Science, Technische Universitat Munchen
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∗)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 136
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 137
Scientific Computing in Computer Science, Technische Universitat Munchen
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 = "[-+][0-9]\.[0-9]+e[-+][0-9]+"
• noch kurzer: rI = "[-+]\d\.\d+e[-+]\d+"
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 138
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 139
Scientific Computing in Computer Science, Technische Universitat Munchen
Regulare Ausdrucke in PythonPattern Syntax• Regulare Ausdrucke immer als ”raw string“• E.g. r’([^.]*.(.*))’. 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 140
Scientific Computing in Computer Science, Technische Universitat Munchen
Maskierungszeichen• Zeichen mit einer speziellen Bedeutung (., *, ...) konnen mit
ihrer eigentlichen Bedeutung durch voranstellen eines\verwendet werden, z.B. \., \*
• Normale Maskierungszeichen (\n, \t, ...) funktionieren wieerwartet, z.B. r’\n+’ entspricht einer 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\W beliebiges nicht-alphanumerisches Zeichen\A entspricht dem Beginn eines strings\Z entspricht dem Ende eines strings
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 141
Scientific Computing in Computer Science, Technische Universitat Munchen
Regular Expressions - erste Python-Beispiele• findall(patt, string) – finde alle nicht uberlappenden
Vorkommen• sub(patt, repl, string) – ersetze Vorkommen durch repl
regstr = r’\d*\.\d*|\d*’
s = "12.3/5/4.4;5.7;6"
findall(regstr ,s)
regstr = r’(\d*\.\d*|\d*)’
sub(regstr , r’\1xxx’, s)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 142
Scientific Computing in Computer Science, Technische Universitat Munchen
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)))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 143
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil VIII
Exceptions
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 144
Scientific Computing in Computer Science, Technische Universitat Munchen
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
Please enter a number: 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 145
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 146
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ...")
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 147
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 148
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 149
Scientific Computing in Computer Science, Technische Universitat Munchen
Build-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]")
except ValueError as e:
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 150
Scientific Computing in Computer Science, Technische Universitat Munchen
Verschachtelte Try-Blocke
def readfloat2 ():
while True:
try:
a = raw_input("Zahl zwischen 0.0 und 1.0: ")
try:
a = float(a)
except ValueError as e:
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]")
except ValueError as e:
print "Fehler: %s" % e
else:
return a
• Exceptions konnen abgefangen und gleich wieder geworfenwerden.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 151
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ():
...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 152
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil IX
2D-Grafiken mit Turtle und TKInter
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 153
Scientific Computing in Computer Science, Technische Universitat Munchen
Zeitslot fuer Abnahme Programmieraufgabe
Zwei mogliche Termine:1) Do., 15.12. von 12:30 - 18:00 Uhr2) Fr., 16.12. von 08:30 - 12:00 Uhr
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 154
Scientific Computing in Computer Science, Technische Universitat Munchen
Visualisierung
R.W. Hamming:”The purpose of computing is insight, not numbers.“
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 155
Scientific Computing in Computer Science, Technische Universitat Munchen
Turtle-Grafiken
• Anschaulicher Beschreibungssprache zum Zeichnen einfacherGrafiken
• Eine Schildkrote lauft uber die Bildflache.• Sie lauft nur gerade Strecken.• Dabei kann sie die Strecke zeichnen (oder auch nicht).• Am Ende angelangt, kann sie sich drehen.• Sehr einfach zu verwenden• Insbesondere fur Fraktale gut geeignet
Einfacher Formalismus• Bewegung um Lange x mit Zeichnung• Bewegung um Lange x ohne Zeichnung• Drehung um Winkel α im Uhrzeigersinn• Drehung um Winkel α gegen den Uhrzeigersinn
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 156
Scientific Computing in Computer Science, Technische Universitat Munchen
Turtle-Grafiken: Modul turtle: erste Befehle• reset(): Zeichenflache erzeugen/zurucksetzen• forward(laenge): laenge Pixel in Blickrichtung laufen• left(winkel): Um winkel Grad nach links drehen• right(winkel): Um winkel Grad nach rechts drehen• goto(x,y): Zu den Koordinaten x,y gehen• penup(): Stift hochnehmen (nicht mehr zeichnen)• pendown(): Stift absenken
import turtle as t
t.reset ()
for i in range (10):
t.forward (50)
t.left (36)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 157
Scientific Computing in Computer Science, Technische Universitat Munchen
Turtle-Grafiken: weitere Befehle• undo() letzten Schritt ruckgangig machen• showturtle(), hideturtle() Schildkrote anzeigen/verstecken• pencolor(color) Farbe auswahlen• circle(radius, extend, steps) Kreis malen• dot(size, color) Punkt malen• speed(n) Malgeschwindigkeit 1-10 wahlen• setworldcoordinates(llx, lly, urx, ury) neue Bezugspunkte• exitonclick() Bild einfrieren bis zum Schließen durch Klick• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 158
Scientific Computing in Computer Science, Technische Universitat Munchen
Chaos: logistische AbbildungIterationsfunktionen / diskrete Abbildungen• xn+1 = Φ(xn): Wert zur Zeit n + 1 ist Funktion des Wertes zur
Zeit n• Verwendung z.B. in der Populationsdynamik• Logistische Abbildung: Φ(x) = rx(1− x)
• z.B. r = 1.5, x0 = 0.1
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 159
Scientific Computing in Computer Science, Technische Universitat Munchen
Zeichnen der Abbildung mit Turtle-Grafik• Datei erstellen (logabb.py)• Turtle einbinden
import turtle
• Die logistische Abbildung:
def log_abb_1(r,x):
return r*x*(1-x)
• Initialisierung
def init(n):
turtle.reset ()
turtle.setworldcoordinates (-1.0,-0.1, n+1, 1.1)
turtle.speed (0)
turtle.hideturtle ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 160
Scientific Computing in Computer Science, Technische Universitat Munchen
• Funktion zum Plotten:
def plot(abb , x_0 , r, numsteps , farbe):
x = x_0
turtle.pencolor(farbe)
turtle.penup ()
turtle.goto(0, x)
turtle.pendown ()
for i in range(numsteps ):
x=abb(r,x)
turtle.goto(i+1,x)
• Plotten fur r = 1.5, x0 = 0.1
from logabb import *
n=100
init(n)
plot(log_abb_1 , 0.1, 1.5, n, "blue")
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 161
Scientific Computing in Computer Science, Technische Universitat Munchen
Ubergang ins Chaos• z.B. r = 3.2, x0 = 0.1
• r = 3.57, r = 3.58, r = 3.59
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 162
Scientific Computing in Computer Science, Technische Universitat Munchen
Auswirkung von Rundungsfehlern• Bisherige Abbildung
def log_abb_1(r,x):
return r*x*(1-x)
• Nun umgeformte Variante (aquivalent!)
def log_abb_2(r,x):
return r*(x-x**2)
• Beide Varianten fur r = 4.0
from logabb import *
n=100
init(n)
plot(log_abb_1 , 0.1, 4.0, n, "blue")
plot(log_abb_2 , 0.1, 4.0, n, "red")
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 163
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 164
Scientific Computing in Computer Science, Technische Universitat Munchen
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• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 165
Scientific Computing in Computer Science, Technische Universitat Munchen
Rekursive Koch-Funktion
import turtle
def kochkurve(laenge , ebene):
if(ebene >0):
kochkurve(laenge /3.0, ebene -1)
turtle.left (60)
kochkurve(laenge /3.0, ebene -1)
turtle.right (120)
kochkurve(laenge /3.0, ebene -1)
turtle.left (60)
kochkurve(laenge /3.0, ebene -1)
else:
turtle.forward(laenge)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 166
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 167
Scientific Computing in Computer Science, Technische Universitat Munchen
Randnotiz: Raumfullende Kurven
• Analoge Erzeugung uber Rekursion• Zusatzlich: Grammatik notig (vgl. re)• Konkrete Nutzung: Parallelisierung, Cache-Nutzung• Beispiele: Hilbert-Kurve, Peano-Kurve
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 168
Scientific Computing in Computer Science, Technische Universitat Munchen
Evaluierung!
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 169
Scientific Computing in Computer Science, Technische Universitat Munchen
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• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 170
Scientific Computing in Computer Science, Technische Universitat Munchen
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
• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 171
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 172
Scientific Computing in Computer Science, Technische Universitat Munchen
Prinzipielle Vorgehensweise• Festlegen, wie die GUI aussehen soll• Zugrundeliegende Funktionalitat implementieren• GUI-Elemente mit Funktionalitat verbinden
Hello world als Klasse
import Tkinter as tk
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 173
Scientific Computing in Computer Science, Technische Universitat Munchen
Anderung der Anordnung
import Tkinter as tk
class MyVis2:
def __init__(self , root):
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 174
Scientific Computing in Computer Science, Technische Universitat Munchen
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 asoziierte Funktion aufgerufen.
class MyVis2:
def __init__(self , root):
...
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 175
Scientific Computing in Computer Science, Technische Universitat Munchen
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:
def __init__(self , root):
self.button1 = tk.Button(root ,
text="Hallo -Knopf", background="blue",
command=self.action1)
self.button2 = tk.Button(root ,
command=self.action2)
...
def action1(self):
...
def action2(self):
...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 176
Scientific Computing in Computer Science, Technische Universitat Munchen
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):
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 177
Scientific Computing in Computer Science, Technische Universitat Munchen
widget Hierarchie
class MyApp:
def __init__(self , root):
# 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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 178
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 179
Scientific Computing in Computer Science, Technische Universitat Munchen
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
• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 180
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 181
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil X
Simulation von Partikelsystemen
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 182
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 183
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 184
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 185
Scientific Computing in Computer Science, Technische Universitat Munchen
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!
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 186
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 187
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 188
Scientific Computing in Computer Science, Technische Universitat Munchen
UML-Diagram (WH)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 189
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 190
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 191
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 192
Scientific Computing in Computer Science, Technische Universitat Munchen
UML-Diagram (WH)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 193
Scientific Computing in Computer Science, Technische Universitat Munchen
Klassen fur Kugeln: Sphere und GravSphere
Konstruktor __init__
• 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 194
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 195
Scientific Computing in Computer Science, Technische Universitat Munchen
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
def calcForce(self , sphere2 ):
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 196
Scientific Computing in Computer Science, Technische Universitat Munchen
UML-Diagram (WH)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 197
Scientific Computing in Computer Science, Technische Universitat Munchen
Schnittstellenklasse zum Steuern der Simulation
• Eine Klasse, die die bisherigen Teile zusammenfuhrt(SimulationInter)
• Dinge, die getan werden mussen:• Initialisierung (Erstellung von Kugeln und Randern)• Erstellung des Integrators• Berechnung der Krafte zwischen den Partikeln• Anwendung der Randbedingungen• 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: GravInter
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 198
Scientific Computing in Computer Science, Technische Universitat Munchen
Notige Module
from spheres import GravSphere , LJSphere
from boundary import Boundary
from spherevis import Spherevis
from Tkinter import Tk, Canvas
from math import sqrt , pi
from random import random
Konstruktor __init__(self, dt)
• Legt Listen fur Kugeln und Rander an (spheres[] undboundaries[])
• Speichert Zeitschrittweite dt
• Erzeugt ein Canvas zur VisualisierungMethode setIntegrator
• Speichert den ubergebenen Integrator ab
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 199
Scientific Computing in Computer Science, Technische Universitat Munchen
Methode calculateForces
• Setzt zunachst die Kraft auf alle Kugeln auf Null• Lauft dann uber alle Kugelpaare und berechnet die Kraft
Methode applyBoundaries
• Pruft fur jede Kugel, ob sie den Rand beruhrt• Falls ja, wird der Effekt des Aufpralls berechnet (Position und
Geschwindigkeit der Kugel angepasst)Methode initVis
• Initialisiert die Visualisierung• Wird durch Subklasse implementiert, da nur da die
Simulationskoordinaten nicht bekannt sindMethode simLoop
• Fuhrt die Simulation endlos durch• In jeden Schritt: Visualisierung, Kraftberechnung, Integration und
Anwendung der Randbedingungen
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 200
Scientific Computing in Computer Science, Technische Universitat Munchen
Implementierung der Klasse SimulationInter
class SimulationInter(object ):
def __init__(self , dt):
self.boundaries = []
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 ()
for sphere1 in self.spheres:
for sphere2 in self.spheres:
if sphere1 != sphere2:
sphere1.calcForce(sphere2)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 201
Scientific Computing in Computer Science, Technische Universitat Munchen
Implementierung der Klasse SimulationInter (cont.)
def applyBoundaries(self):
for b in self.boundaries:
for p in self.spheres:
b.boundaryAction(p)
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)
self.applyBoundaries ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 202
Scientific Computing in Computer Science, Technische Universitat Munchen
UML-Diagram (WH)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 203
Scientific Computing in Computer Science, Technische Universitat Munchen
Subklasse AstroInter
Konstruktor __init__
• 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))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 204
Scientific Computing in Computer Science, Technische Universitat Munchen
Restliche Planeten
self.spheres.append(GravSphere (( 57.9e9, 0.),
(0. ,47872.) , c*10.e9, 3.302 e23)) # Merkur
self.spheres.append(GravSphere (( 108.2e9, 0.),
(0. ,35021.) , c*10.e9, 4.8685 e24)) # Venus
self.spheres.append(GravSphere (( 149.6e9, 0.),
(0. ,29786.) , c*10.e9, 5.9736e24 , "blue")) # Erde
self.spheres.append(GravSphere (( 227.9e9, 0.),
(0. ,24131.) , c*10.e9, 6.4185 e23)) # Mars
self.spheres.append(GravSphere (( 778.3e9, 0.),
(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
self.spheres.append(GravSphere ((2.875e12 , 0.),
(0., 6835.) , c*10.e9, 8.6832 e25)) # Uranus
self.spheres.append(GravSphere ((4.504e12 , 0.),
(0., 5478.) , c*10.e9, 1.0243 e26)) # Neptun
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 205
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 206
Scientific Computing in Computer Science, Technische Universitat Munchen
Screenshot der Astro-Simulation
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 207
Scientific Computing in Computer Science, Technische Universitat Munchen
Mogliche Erweiterung?
Welche anderen physikalischen Phanomene lassen sich ahnlichbeschreiben?
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 208
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 209
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 210
Scientific Computing in Computer Science, Technische Universitat Munchen
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))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 211
Scientific Computing in Computer Science, Technische Universitat Munchen
Erweiterung auf MD-Simulation
Euler
+ integrate
Integrator- dt: float
+ integrate
Spherevis- dt: float
- pos2pixel+ drawSpheres
Boundary- dim, pos, side
- distance- particleTouches+ boundaryAction
SimulationInter+ spheres+ boundaries+ canvas, dt
- calculateForces- applyBoundaries+ setIntegrator+ initVis+ simLoop
GravSphere
+ calcForce
Sphere+ rad, mass, pos+ vel, force, color- dist+ resetForce+ calcForce
LJSphere
+ calcForce
AstroInter
+ initVis
MDInter
+ initVis
n
n 1
1
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 212
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ):
def calcForce(self , sphere2 ):
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]
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 213
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 214
Scientific Computing in Computer Science, Technische Universitat Munchen
Klasse boundary
Konstruktor __init__
• 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 215
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 216
Scientific Computing in Computer Science, Technische Universitat Munchen
Subklasse MDInter
Konstruktor __init__
• 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 217
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 218
Scientific Computing in Computer Science, Technische Universitat Munchen
Einbinden der neuen Klassennotige Module
from optparse import OptionParser
from integrator import Euler
from simulationinterface import AstroInter , MDInter
Option-Parser
# Optionen festlegen
p = OptionParser ()
p.add_option("-t", "--simType", dest="simType",
help="’Astro ’ or ’MD’ simulation")
p.add_option("-n", type="int", dest="numParts",
help="Number of particles", default =30)
p.add_option("--rho", type="float", dest="rho",
help="Density", default =0.1)
(options , args) = p.parse_args ()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 219
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 220
Scientific Computing in Computer Science, Technische Universitat Munchen
Screenshot der MD-Simulation
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 221
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil XI
Wissenschaftliches Rechnen in Python
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 222
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 223
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 224
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 225
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 226
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 227
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
• Ausnahme: object_ speichert Zeiger auf Objekte
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 228
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 229
Scientific Computing in Computer Science, Technische Universitat Munchen
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))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 230
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 231
Scientific Computing in Computer Science, Technische Universitat Munchen
Lesen und Schreiben von Array-Elementen
• Die von Listen bekannte Indizierung kann verwendet werden• Dimensionen werden mit Komma separiert
a = arange (20); a.resize ((4 ,5))
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 232
Scientific Computing in Computer Science, Technische Universitat Munchen
Uber Elemente von Matrizen iterieren
a = arange (20); a.resize ((4 ,5))
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 233
Scientific Computing in Computer Science, Technische Universitat Munchen
Arbeiten mit Arrays• Schnelle built-in Methoden fur Arrays
a = arange (12); a.resize ((3 ,4))
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 234
Scientific Computing in Computer Science, Technische Universitat Munchen
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.I # Inverse
m.A # Als 2d Array
m.H # Konjugiert Transponierte
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 235
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 236
Scientific Computing in Computer Science, Technische Universitat Munchen
Untermodule (2)
Modul numpy.linalg• Die wichtigsten Werkzeuge fur lineare Algebra
>>> from numpy.linalg import *
>>> a = arange (12); a.resize ((3 ,4))
>>> 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 . . .
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 237
Scientific Computing in Computer Science, Technische Universitat Munchen
Matplotlib
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 238
Scientific Computing in Computer Science, Technische Universitat Munchen
MatplotlibWozu 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 aufVersion-Durcheinander• Numpy, scipy, ipython und matplotlib werden oft gleichzeitig
benutzt• Die Pakete hangen voneinander ab (matplotlib benutzt z.B.
numpy-Arrays• Pylab ist ein (inoffizielles) Paket, das alle vier enthalt• Je nach Betriebssystem(sversion) mussen evtl. unterschiedliche
Pakete installiert werden (D.h. Der Modulname im import-Befehlkann auch unterschiedlich sein)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 239
Scientific Computing in Computer Science, Technische Universitat Munchen
Moglichkeiten zum importierenpython/ipython interactive
> ipython
import scipy , matplotlib.pylab
x = scipy.randn (10000)
matplotlib.pylab.hist(x, 100)
> ipython
import numpy.random , matplotlib.pylab
x = numpy.random.randn (10000)
matplotlib.pylab.hist(x, 100)
ipython mit pylab-Modus
> ipython -pylab
x = randn (10000)
hist(x, 100)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 240
Scientific Computing in Computer Science, Technische Universitat Munchen
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()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 241
Scientific Computing in Computer Science, Technische Universitat Munchen
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.’)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 242
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ’)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 243
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel - Histogramm
# entweder ipython -pylab oder folgende imports:
#from numpy import *
#from matplotlib.mlab import *
#from matplotlib.pyplot import *
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()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 244
Scientific Computing in Computer Science, Technische Universitat Munchen
SciPy
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 245
Scientific Computing in Computer Science, Technische Universitat Munchen
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 (teils aufbauend auf NumPy und SciPy)• Hier daher nur eine kurze Ubersicht• Weitere Informationen auf www.scipy.org
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 246
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 247
Scientific Computing in Computer Science, Technische Universitat Munchen
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-Function• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 248
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Interpolation - Linear
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
x = np.arange (0,10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y)
xnew = np.arange (0 ,9,0.1)
plt.plot(x,y,’o’,xnew ,f(xnew),’-’)
plt.title(’Lineare Interpolation ’)
plt.show()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 249
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Interpolation - Kubische Splines
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)
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()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 250
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil XII
Datenstrukturen: Baume, Stacks undQueues
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 251
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 252
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 253
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
hanoi(3,A,B,H)
hanoi(2,A,H,B)
hanoi(1,A,B,H)
A -> B
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 254
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
hanoi(3,A,B,H)
hanoi(2,A,H,B)
hanoi(1,A,B,H)
A -> B
A -> H
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 255
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 256
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 257
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 258
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 259
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 260
Scientific Computing in Computer Science, Technische Universitat Munchen
A B Hilfe
def hanoi(n, a, b, hilfe):
if(n > 0):
hanoi(n-1, a, hilfe , b)
b.append(a.pop())
hanoi(n-1, hilfe , b, a)
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 261
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 262
Scientific Computing in Computer Science, Technische Universitat Munchen
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)• ...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 263
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: Routensuche
Abbildung: Dijkstra bzw. A*; Quelle: Bungartz, Zimmer, Buchholz, Pfluger:Modellbildung und Simulation - Eine anwendungsorientierte Einfuhrung
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 264
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 265
Scientific Computing in Computer Science, Technische Universitat Munchen
(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
IN8008, Wintersemester 2011/2012 266
Scientific Computing in Computer Science, Technische Universitat Munchen
Tiefensuche in (unsortierten) Baumen/Graphen• Erster Nachbarknoten, dann dessen erster 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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 267
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 268
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel Breiten- und Tiefensuche
1
2 3
4 75 6
8 9
1
7 2
9 38 4
6 5
Breitensuche Tiefensuche
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 269
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 270
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 271
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 272
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 273
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 274
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 275
Scientific Computing in Computer Science, Technische Universitat Munchen
Loschen eines Elements• Vor dem Loschen sind alle Balancewerte ∈ −1,0,1.• Zunachst wir das zu loschende Element 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 genau einen Sohn s• 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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 276
Scientific Computing in Computer Science, Technische Universitat Munchen
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’
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 277
Scientific Computing in Computer Science, Technische Universitat Munchen
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!
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 278
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ’ ’
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 279
Scientific Computing in Computer Science, Technische Universitat Munchen
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))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 280
Scientific Computing in Computer Science, Technische Universitat Munchen
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’])
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 281
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 282
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 283
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil XIII
Simulation mit PDEs:Warmeleitungsgleichung
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 284
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 285
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 286
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 287
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 288
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel: WarmeleitungKernproblem der Thermodynamik• Warme beeinflusse die Oberflache eines Objekts• Ziel: Ausbreitung bzw. Verteilung 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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 289
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 290
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 291
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 292
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 293
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 294
Scientific Computing in Computer Science, Technische Universitat Munchen
Bilanzgleichung• Anderung der Warme entspricht Zu- und Abflussen
dQdt
=
∫V
cρ∂T∂t
(~x ; t) d~x =
∫∂V
k∇T (~x ; t)~n−→dS
• Gaußscher Integralsatz∫V
cρ∂T∂t
(~x ; t) d~x =
∫V
k∆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) ,
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 295
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 296
Scientific Computing in Computer Science, Technische Universitat Munchen
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
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 297
Scientific Computing in Computer Science, Technische Universitat Munchen
Adaptivitat• Feineres Gitter⇒ teuerer und genauer• Oft ist die hohe Genauigkeit nur an wenigen Stellen notig
• Strukturmechanik: ”Sollbruchstellen“• Stromungsmechanik: Wirbel• Warmeleitung: dunne Bauteile
Ansatze• Finite Differenzen• Finite Volumen• Finite Elemente
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 298
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 299
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 300
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 301
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 302
Scientific Computing in Computer Science, Technische Universitat Munchen
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
=
−T0
0...0
−Tn+1
.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 303
Scientific Computing in Computer Science, Technische Universitat Munchen
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 .
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 304
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 305
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 306
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil XIV
Losung linearer Gleichungssysteme
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 307
Scientific Computing in Computer Science, Technische Universitat Munchen
Gauß-Algorithmus
• Zwei Schritte: Vorwartselimination und RuckwartssubstitutionVorwartselimination• Erzeugen einer Stufenform• Zeilen durfen mit Konstanten multipliziert werden• Zeilen durfen addiert werden werden a11 a12 a13
a21 a22 a23a31 a32 a33
· x1
x2x3
=
b1b2b3
- 1. Zeile · a21a11
- 1. Zeile · a31a11 a11 a12 a13
0 a22 a230 a32 a33
· x1
x2x3
=
b1
b2
b3
- 2. Zeile · a32
a22 a11 a12 a130 a22 a230 0 a33
· x1
x2x3
=
b1
b2
b3
- 2. Zeile · a32
a22
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 308
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 309
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel
2 3 14 7 32 4 4
· x1
x2x3
=
122
2 3 1
0 1 10 1 3
· x1
x2x3
=
101
2 3 1
0 1 10 0 2
· x1
x2x3
=
101
x3 = 1/2 = 0.5x2 = (0− 1 · 0.5)/1 = −0.5x1 = (1− 3 · (−0.5)− 1 · 0.5)/2 = 1.0
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 310
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 311
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 312
Scientific Computing in Computer Science, Technische Universitat Munchen
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])
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 313
Scientific Computing in Computer Science, Technische Universitat Munchen
Losung der Warmeleitungsgleichung 1D
• Stab mit neun Diskretisierungspunkten• Links beheizt: T = 1, rechter Rand: T = 0
>>> def getwaermelgs1D ():
>>> a = array([[-2, 1, 0, 0, 0, 0, 0, 0, 0],
>>> [ 1,-2, 1, 0, 0, 0, 0, 0, 0],
...
>>> [ 0, 0, 0, 0, 0, 0, 0, 1,-2]],
>>> dtype = float)
>>> b = array([-1, 0, 0, 0, 0, 0, 0, 0, 0],
>>> dtype = float)
>>> return(a,b)
>>>
>>> (a,b) = getwaermelgs1D ()
>>> loese(a,b)
array ([0.9 , 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1])
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 314
Scientific Computing in Computer Science, Technische Universitat Munchen
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])
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 315
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 316
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 317
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 318
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 319
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 320
Scientific Computing in Computer Science, Technische Universitat Munchen
Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 321
Scientific Computing in Computer Science, Technische Universitat Munchen
Jacobi-Verfahren (aus PDE-Sicht)• Ist ein Relaxationsverfahren• An jedem Punkt wird 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 322
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 323
Scientific Computing in Computer Science, Technische Universitat Munchen
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
IN8008, Wintersemester 2011/2012 324
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 325
Scientific Computing in Computer Science, Technische Universitat Munchen
ImplementierungMaß 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 *
import matplotlib.pyplot as plt
from math import sqrt
class Waermesimulation:
...
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 326
Scientific Computing in Computer Science, Technische Universitat Munchen
ImplementierungKonstruktor• 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 327
Scientific Computing in Computer Science, Technische Universitat Munchen
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 i in range(1, self.zeilen +1):
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 ))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 328
Scientific Computing in Computer Science, Technische Universitat Munchen
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()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 329
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 330
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 331
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T7
T4
T6
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 332
Scientific Computing in Computer Science, Technische Universitat Munchen
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
for i in range(1, self.zeilen +1):
for j in range(1, self.spalten +1):
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 ))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 333
Scientific Computing in Computer Science, Technische Universitat Munchen
Teil XV
Mehrgitterverfahren zur Losung linearerGleichungssysteme
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 334
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 335
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 336
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uu uu 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 337
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uu u uu 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 338
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uu uu
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 339
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uuu uu
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 340
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uuu u uu
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 341
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uuu uu 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 342
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uuu 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 343
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uuu u 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 344
Scientific Computing in Computer Science, Technische Universitat Munchen
u uuu uuuu uuu 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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 345
Scientific Computing in Computer Science, Technische Universitat Munchen
Konvergenz: Fourier-Analyse
Szenario• Stationare Warmeleitung bei quadratischer Flache• Randbedingung uberall Null• Losung: Null auf der gesamten Flache• Initialer x-Vektor ist damit gleich dem Fehler e• Veranderung des Fehlers in einer Iteration wird untersucht• Initialer Fehler als Kombination von Sinus-Kurven
unterschiedlicher Frequenz• 1D: x =
∑m am · (sin(mπhi))i=1..N
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 346
Scientific Computing in Computer Science, Technische Universitat Munchen
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!
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 347
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 348
Scientific Computing in Computer Science, Technische Universitat Munchen
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 Gitterpunkte notig• ⇒ Anzahl an Iterationen wachst mit 1/h2
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 349
Scientific Computing in Computer Science, Technische Universitat Munchen
Frequenzabhangige Fehlerreduktion
Jacobi-Verfahren
πh π2hπNh
cos
π
gedampftes Jacobi-Verfahren
πh π2hπNh π
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 350
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 351
Scientific Computing in Computer Science, Technische Universitat Munchen
Mehrgitter: Motivation
• Optimales LGS-Verfahren: Rechenaufwand wachst linear mit derAnzahl ab Diskretisierungspunkten
• Besser geht es nicht (jeder Punkt muss mindestens einmalbetrachtet werden)
• Kosten einer Iteration bei Jacobi/Gauss-Seidel sind O(n)
• Die Anzahl an Iterationen wachst mit der Anzahl anDiskretisierungspunkten
• 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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 352
Scientific Computing in Computer Science, Technische Universitat Munchen
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.
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 353
Scientific Computing in Computer Science, Technische Universitat Munchen
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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 354
Scientific Computing in Computer Science, Technische Universitat Munchen
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 Punkte wird
verwendet, die dazwischen verworfen• Teilweise bessere Ergebnisse, wenn zwischen drei
benachbarten Punkten gemittelt wird
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 355
Scientific Computing in Computer Science, Technische Universitat Munchen
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
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 356
Scientific Computing in Computer Science, Technische Universitat Munchen
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
self.glaetter(2, level)
# Restriktion: nicht notig (Injektion)
# -
# Loesung auf grobem Gitter
self.mehrgitter_rekursiv(level +1)
# Prolongation
self.prolongation(level)
# Nachglaetten / oberster Level
self.glaetter(2, level)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 357
Scientific Computing in Computer Science, Technische Universitat Munchen
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 ])
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 358
Scientific Computing in Computer Science, Technische Universitat Munchen
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
for i in range(1, self.zeilen +1):
for j in range(1, self.spalten +1):
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 ))
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 359
Scientific Computing in Computer Science, Technische Universitat Munchen
Beispiel
>>> from waermesim import *
>>> lgs = Waermesimulation (50. ,30. ,1. ,0.6 ,0. ,1. ,0.)
>>> lgs.mehrgitter (0.001)
>>> lgs.plot()
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 360
Scientific Computing in Computer Science, Technische Universitat Munchen
Kosten des Mehrgitter-Verfahrens
• 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)
T. Neckel: Einfuhrung in die wissenschaftliche Programmierung
IN8008, Wintersemester 2011/2012 361