Wissenschaftliches Rechnen mit Python

Zweistündiger Workshop am Linux Infotag 2013, 23. März 2013

Hochschule Augsburg

Hubert Högl, <Hubert.Hoegl@hs-augsburg.de>

URL: http://elk.informatik.fh-augsburg.de/pub/LIT/2013/scipy-workshop/index.html

Inhalt

1 Was muss installiert werden?
2 Programme laufen lassen
3 IPython
4 Ein Python Programm zur Erinnerung
5 Einfache Rechnungen mit der Standardbibliothek
6 NumPy
7 Matplotlib
8 Mayavi
9 SciPy
10 SymPy
11 Aus folgenden Texten wurde zitiert
12 Mehr zum Thema

1 Was muss installiert werden?

Linux Anwender installieren folgende Pakete: vim, ipython, ipython-qtconsole, ipython-notebook, python-numpy, python-scipy, python-matplotlib, python-sympy, spyder (oder python-spyderlib), mayavi2.

Das ist nur eine kleine Auswahl aller vorhandenen Pakete unter Linux. Im wesentlichen lassen sich alle Pakete, die bei Python(x, y) schon gepackt sind, auch unter Linux installieren. Einige Hinweise enthält auch der Text "Py4Science" (http://fperez.org/py4science/starter_kit.html) von Fernando Perez, dem IPython Autor.

Noch ein Hinweis: Die Pakete ipython-qtconsole und ipython-notebook gibt es bei Debian (Squeeze) nicht. Bei Ubuntu Varianten (Xubuntu, Mint) sind sie dabei.
Windows Anwender installieren das Python(x,y) Paket, siehe http://code.google.com/p/pythonxy.

2 Programme laufen lassen

Alternativen

Interaktiver Python Prompt
```
python [RETURN]

Python 2.7.3 (default, Sep 26 2012, 21:53:58)
Type "help", "copyright", "credits" or "license" for more information.
>>> f = 4.2
>>> f
4.2
```
Beachte:
- Navigation mit Emacs Tasten: C-P, C-N, C-A, C-E, C-B, C-F, C-K, C-Y
- Vorsicht: Eingaben gehen verloren.
- IPython bietet viel mehr Möglichkeiten.
Quelltextdatei vom Interpreter ausführen lassen.
```
python programm.py [RETURN]
```

3 IPython

Interactive Python

hhoegl@aspire ~/TempGit/python-arbeitsbuch $ ipython
Python 2.7.3 (default, Sep 26 2012, 21:53:58)
Type "copyright", "credits" or "license" for more information.

IPython 0.13.1.rc2 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]:

IPython kann man auf verschiedene Weise starten:

Mit Zeichen-Schnittstelle: ipython [RETURN]

Mit GUI-Schnittstelle: ipython qtconsole [RETURN]

Mit Web-Schnittstelle: ipython notebook [RETURN]

Aufgaben

IPython starten und mit ? den einführenden Text aufrufen.
Mit quickref die vorhandenen speziellen Kommandos ("magic commands") ausgeben.

Wichtige Eigenschaften der IPython Shell sind:
- Hilfe (name?, name??)
- History (_i, _ii, _iii,_ih[n])
- Schnelle Navigation in Verzeichnissen (%dhist, %dirs, %cd, %ls)
- Programme laufen lassen (%run)
- Debuggen (%run -d)
- Zeiten messen (%time, %timeit)
- Logging (%logon, %logoff, %logstart)
- Editieren (%edit)
- viele mehr...

4 Ein Python Programm zur Erinnerung

# -*- coding: utf-8 -*-

'''
Erste Übung (eins.py)
'''

TEILNEHMER =  10
for i in range(TEILNEHMER):
    print "Hallo Workshopteilnehmer %d!" % i

def umfang(radius):
    '''Diese Funktion berechnet den Umfang
    eines Kreises.
    radius -- Der Radius des Kreises
    '''
    import math
    return 2 * radius * math.pi

r = 2.0
u = umfang(r)

print
print "Bei einem Radius von %1.2f ist der Umfang %1.2f." % (r, u)

Code: workshop/eins.py

Python 2
Kodierung
Docstrings
Variablenzuweisung (Konvention: Gross = Konstante)
Blöcke 4 Leerzeichen einrücken!
Datentypen:
- None
- Boolean: True, False
- Integer: 42
- Gleitkommazahl: 7.34
- Komplexe Zahl: 3+6j
- String: "Hallo" ("Sequenztyp")
- Liste: [1, 10, 100, 1000] ("Sequenztyp")
- Tupel: (1, 10, 100, 1000) ("Sequenztyp")
- Dictionary: {'rot' : 0xff0000, 'gruen' : 0x00ff00, 'blau' : 0x0000ff}
Range gibt Sequenztyp zurück
Formatstring Anweisungen wie in C
Funktion, math aus Python Stdlib ("batteries included")

5 Einfache Rechnungen mit der Standardbibliothek

Mathematik: math, cmath ("complex math")
Keine Arrays, keine Matrizen!

Komplexe Zahlen

import math

print math.sqrt(2)
# print math.sqrt(-2)
print

import cmath
print cmath.sqrt(2)
print cmath.sqrt(-2)
print

import numpy.lib.scimath as sm
print sm.sqrt(2)
print sm.sqrt(-2)

Code: workshop/complex.py

Aufgaben

Kurz den Funktionsumfang der Module math und cmath überfliegen:

http://docs.python.org/2.7/library/index.html
Prüfen Sie, ob folgende Formel stimmt:

\begin{equation*} e^{iq} = cos(q) + i*sin(q) \end{equation*}
Fluglinie eines Balles ausgeben. [LANGTANGEN]

\begin{equation*} f(x) = x \, tan ( \theta ) - \frac{1}{2 v^2_o} \frac{gx^2}{cos^2(\theta)} + y_o \end{equation*}

workshop/flugbahn.py
Der Widerstand R eines BetaTherm Temperatursensors gehorcht der Steinhart-Hart Gleichung. Diese lautet:

\begin{equation*} a = (A - (1/t)) / (2 * C) \end{equation*}

\begin{equation*} b = \sqrt{ \frac{((A - (1/t))/C)^{2}}{4} + \frac{(B/C)^3}{27}} \end{equation*}

\begin{equation*} c = \sqrt[3]{ -(a + b) } + \sqrt[3]{ -(a - b) } \end{equation*}

\begin{equation*} R = e^{c} \end{equation*}

Die Temperatur t muss in Kelvin angegeben werden. Man kommt von den Einheiten Grad auf Kelvin, wenn man 273.13 dazuzählt.

Die Konstanten A, B und C lauten: A = 1.471388E-3, B = 2.376138E-4, C = 1.051058E-7.

Ihre Aufgabe ist nun, den Widerstandsverlauf zwischen -20 und 140 Grad Celsius zu berechnen.

Hinweis: Es gibt keine Kubikwurzel in Python. Verwenden Sie deshalb die pow(x, y) Funktion mit gebrochenem Exponenten.

Lit.: http://www.beta.dk/betathermkatalog/teoridel.pdf

6 NumPy

Grosse multi-dimensionale Arrays und Matrizen, mit Bibliothek für high-level Funktionen.

Grobe Struktur von NumPy:

Array objects
- N-dim array (ndarray)
- scalars
- Data type objects (dtype)
- Indexing
- ...
Universal functions (ufunc)
- Broadcasting
Routines
- Array creation and manipulation, indexing
- input/output
- numpy.fft
- numpy.linalg
- numpy.random
- Sorting, searching, and counting
- Statistics
- Mathematical functions
- Functional programming
- Polynomials
- Financial functions
- Matrix library (numpy.matlib)
- ...

Historische Entwicklung: Numeric, numarray, numpy.

ndarray: n-dimensional Array (kurz: array)
Axes: Dimensionen
Rank: Anzahl der Achsen
Shape: Grösse des Arrays in jeder Dimension (Tuple).
SciPy: "Scientific Python" (SciPy) baut auf NumPy auf.

Alle elementaren mathematischen Funktionen die man aus den Standardbibliotheksfunktionen math und cmath kennt, sind auch vorhanden im Modul scimath, das man so importieren kann:

import numpy.lib.scimath as sm

Diese Funktionen geben je nach Resultat-Typ entweder float- oder complexwertige Resultate zurück.

Am Anfang eines Programmes sollte man sich immer einen kurzen Alias der benötigten wissenschaftlichen Module beschaffen. Folgendes wird häufig verwendet:

import numpy as np
import scipy as sp
import matplotlib as mpl
import matplotlib.pyplot as plt

Aufgabe Vollziehen Sie die folgenden Array und Matrizen Eingaben auf Ihrem Python oder IPython Prompt mit geänderten Werten nach.

Arrays

>>> a = np.array((1, 2, 3))
>>> a
array([1, 2, 3])

>>> b = np.array([1, 2, 3])
>>> b

>>> c = np.array([1, 2, 3], dtype=np.complex)
>>> c
array([ 1.+0.j,  2.+0.j,  3.+0.j])

>>> a1 = np.arange(10)
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

>>> a2 = np.arange(10.0)
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.]

>>> a3 = np.arange(1.0, 2.0, 0.25)
>>> a3
array([ 1.  ,  1.25,  1.5 ,  1.75])

>>> d = np.linspace(0, 2, 9)
>>> d
array([ 0.,  0.25,  0.5,  0.75,  1., 1.25, 1.5, 1.75, 2.  ])

>>> a = np.arange(12)
>>> b = a.reshape(3, 4)
>>> b
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

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

>>> a = np.arange(4)
>>> a.resize(2,3)
>>> a
array([[0, 1, 2],    # aufgefuellt mit Nullen!
       [3, 0, 0]])

Array mit Nullen

>>> np.zeros([4, 4])
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])

Array mit Einsen

>>> np.ones([4, 4])
array([[ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.]])

Array mit Einsen in Diagonale

>>> np.eye(4)
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

Index Operator

>>> m = np.matrix([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]])

>>> m.shape
(2, 5)

>>> m[0]
matrix([[1, 2, 3, 4, 5]])

>>> m[1]
matrix([[6, 7, 8, 9, 10]])

>>> m[1, 3]
9

>>> m[0:2, 3]   # equivalent: m[:, 3]
matrix([[4],
        [9]])

>>> m[1, ...]
x([[ 6,  7,  8,  9, 10]])

Rechnen mit Arrays

>>> v = np.array([0, 1, 2, 3, 4])
>>> w = np.array([0, 1, 2, 3, 4])

>>> 5 + v
array([5, 6, 7, 8, 9])

>>> 2 * v
array([0, 2, 4, 6, 8])

>>> v + w
array([0, 2, 4, 6, 8])

>>> v * w
array([ 0,  1,  4,  9, 16])

>>> np.inner(v, w)
30

>>> np.outer(v, w)
array([[ 0,  0,  0,  0,  0],
       [ 0,  1,  2,  3,  4],
       [ 0,  2,  4,  6,  8],
       [ 0,  3,  6,  9, 12],
       [ 0,  4,  8, 12, 16]])

Matrizen aus Arrays

\begin{equation*} m1 = \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ \end{pmatrix} \end{equation*}

\begin{equation*} m2 = \begin{pmatrix} 2 & 1 \\ 4 & 3 \\ \end{pmatrix} \end{equation*}

>>> m1 = np.array([[1, 2], [3, 4]])
>>> m2 = np.array([[2, 1], [4, 3]])

>>> 5 + m1
array([[6, 7],
       [8, 9]])

>>> 2 * m1
array([[2, 4],
       [6, 8]])

>>> m1 + m2        # elementweise Addition
array([[3, 3],
       [7, 7]])

>>> m1 * m2        # elementweise Multiplikation
array([[ 2,  2],
       [12, 12]])

>>> np.dot(m1, m2)
array([[10,  7],   # 'richtige' Matrixmultiplikation
       [22, 15]])

Eingebaute Matrizen

>>> m1 = np.matrix([[1, 2], [3, 4]])
>>> m2 = np.matrix([[2, 1], [4, 3]])

>>> 5 + m1
matrix([[6, 7],
        [8, 9]])

>>> 2 * m1
matrix([[2, 4],
        [6, 8]])

>>> m1 + m2       # elementweise Addition
matrix([[3, 3],
        [7, 7]])

>>> m1 * m2       # richtige Matrixmultiplikation
matrix([[10,  7],
        [22, 15]])

>>> m = np.matrix([[1, 2], [3, 4]])
>>> m.I
matrix([[-2. ,  1. ],
        [ 1.5, -0.5]])

Kurzschreibweise

>>> m1 = np.mat('1, 2; 3, 4')
>>> m1
matrix([[1, 2],
        [3, 4]])

Kreuzprodukt (2-dim, 3-dim)

\begin{equation*} a = \begin{pmatrix} 0, & 1, & 2 \\ \end{pmatrix} \end{equation*}

\begin{equation*} b = \begin{pmatrix} 1, & 3, & 2\\ \end{pmatrix} \end{equation*}

\begin{equation*} \begin{pmatrix} a_2b_3 - a_3b_2 \\ a_3b_1 - a_1b_3 \\ a_1b_2 - a_2b_1 \\ \end{pmatrix} = \vec{a} \times \vec{b} \end{equation*}

>>> v = np.array([0, 1, 2,])
>>> w = np.array([1, 3, 2,])
>>> np.cross(v, w)
array([-4,  2, -1])

Matrix mal Vektor

>>> m = np.array([[1, 2], [3, 4]])
>>> v = np.array([1, 2])
>>> m * v
array([[1, 4],
       [3, 8]])

>>> m = np.matrix([[1, 2], [3, 4]])
>>> v = np.matrix([1, 2])
>>> m * v.T
matrix([[ 5],
        [11]])

Summe, Mittelwert

>>> a = np.arange(0, 101, 1)
>>> a.sum()
5050
>>> a.mean()
50.0

Weitere Operationen der linearen Algebra

trace()
linalg.det()
linalg.norm()
linalg.eigvals()
linalg.solve()
linalg.inv()

Aufgabe: Bestimmen Sie die Determinate der folgenden Matrizen:

\begin{equation*} A = \begin{pmatrix} 2 & -3 \\ 5 & 1 \\ \end{pmatrix} \qquad B = \begin{pmatrix} 1 & 3 & 5 \\ -1 & 2 & 0 \\ 4 & 2 & -3 \\ \end{pmatrix} \end{equation*}

Vektorisierte Funktionen

Elementweise Anwendung auf Array

a = np.linspace(0, 2, 10)
>>> def f(x):
...     return x+1
...
>>> f(a)

Aufgaben Vektorisieren Sie die im obigen Kapitel "Einfache Rechnungen mit der Standardbibliothek" angegebenen Aufgaben (Flugbahn, etc.).

Zufallszahlen

>>> np.random.rand(4, 2)
array([[ 0.53002771,  0.95503893],
       [ 0.416452  ,  0.40608664],
       [ 0.01201623,  0.37152102],
       [ 0.96669925,  0.68950739]])

Literatur

NumPy User Guide

http://docs.scipy.org/doc/numpy/user/index.html

7 Matplotlib

Matplotlib ist eine 2D Plot-Bibliothek. Die Homepage ist hier:

http://matplotlib.org

Es gibt verschiedene "Backends" für Tkinter (TkAgg), Qt4 (QT4Agg), wxWidgets (WXAgg) und Gtk (GTKAgg).

PyLab

Bestimmte Sicht auf NumPy, SciPy und Matplotlib, die es Anfängern besonders einfach machen soll, damit zu arbeiten.

Die Homepage http://www.scipy.org/PyLab sagt:

The point is that PyLab is a compelling, integrated, usable and superior alternative to MATLAB.

In Python importiert man:

import pylab

oder

from pylab import *

Mit help(pylab) bekommt man einen Hilfetext (Docstring).

Man kann ipython mit --pylab starten. Damit wird das TkAgg Backend zur Ausgabe verwendet. Mit --pylab=qt wird QT4Agg zur Ausgabe verwendet.

Im Wesentlichen sind bei PyLab die Namensräume von matplotlib.pyplot (http://matplotlib.org/api/pyplot_api.html) und numpy kombiniert, was bei der interaktive Arbeit Zeit spart. Zum Programmieren sollten aber die Namensräume getrennt werden:

import numpy as np
import matplotlib.pyplot as plt

Einfaches Beispiel:

In [1]: y = [1, 2, 3]
In [2]: plot(y)        # x-Punkte sind index Werte von y

Plots

Beispiel 1
- workshop/plot1.py, workshop/plot1.png
Beispiel 2
- workshop/barchart.py
Beispiel 3

Aufgabe: Suche den Schnittpunkt der zwei Geraden

\begin{equation*} \; \frac{1}{2}x - y = -1 \end{equation*}

\begin{equation*} - \frac{1}{4} - y = -5 \end{equation*}

\begin{equation*} \Rightarrow K = \begin{pmatrix} 0.5 & -1 \\ -0.25 & -1 \\ \end{pmatrix} \qquad B = \begin{pmatrix} -1 \\ -5 \\ \end{pmatrix} \end{equation*}

workshop/geradenschnitt.py, workshop/geradenschnitt.png

Aufgabe Nachdem Sie die oben angegebene Aufgaben (Flugbahn, etc.) zunächst mit der Standardbibliothek gelöst, dann mit NumPy vektorisiert haben, können Sie die Ergebnisse nun mit der Matplotlib grafisch ausgeben.

Aufgabe Sehen Sie sich das "Matplotlib Tutorial" und die "Example List" an und experimentieren Sie mit Plots die Sie schon lange mal machen wollten.

Literatur

Plotting commands summary

http://matplotlib.org/api/pyplot_summary.html
Pyplot API

http://matplotlib.org/api/pyplot_api.html
Matplotlib Tutorial

http://www.loria.fr/~rougier/teaching/matplotlib
Matplotlib example list

http://matplotlib.org/examples/index.html

8 Mayavi

Mayavi ist eine 3D Plot-Software. Die Homepage ist hier:

http://code.enthought.com/projects/mayavi/

Beispiel:

Entnommen aus [MAYAVI1]:

IPython starten mit
```
ipython --pylab=qt
```
oder
```
ipython qtconsole
```
Programm workshop/mlabdemo.py laden [MAYAVI1].

9 SciPy

Subpackage	Description
cluster	Clustering algorithms
constants	Physical and mathematical constants
fftpack	Fast Fourier Transform routines
integrate	Integration and ordinary differential equation solvers
interpolate	Interpolation and smoothing splines
io	Input and Output
linalg	Linear algebra
maxentropy	Maximum entropy methods
ndimage	N-dimensional image processing
odr	Orthogonal distance regression
optimize	Optimization and root-finding routines
signal	Signal processing
sparse	Sparse matrices and associated routines
spatial	Spatial data structures and algorithms
special	Special functions
stats	Statistical distributions and functions
weave	C/C++ integration

Beispiele

optimize.fsolve

[LANGTANGEN] workshop/algebraic.py, workshop/algebraic.png

Aufgabe Suchen Sie die Nullstelle in folgender Kurve:

workshop/algebraic2.py, workshop/algebraic2.png

Aufgabe Holen Sie sich Beispielprogramme aus einem wissenschaftichen Bereich, der Sie interessiert und vollziehen Sie diese nach. Gute Quellen sind das SciPy Kochbuch http://www.scipy.org/Cookbook und die Python Scientific Lecture Notes http://scipy-lectures.github.com.

10 SymPy

SymPy ist eine Python Bibliothek für Symbolische Mathematik (http://sympy.org).

Ein Beispiel:

>>> from sympy import *
>>> x = Symbol('x')
>>> diff(sin(x), x)
cos(x)
>>> y = Symbol('y')
>>> simplify((x+x*y)/x)
y + 1

Aufgaben

Bilden Sie die erste Ableitung der weiter oben im Abschnitt "SciPy" angegebenen algebraischen Gleichungen, bei denen wir die Nullstellen gesucht haben.

Literatur