Physik der sozio-ökonomischen Systeme mit dem Computer¶

(Physics of Socio-Economic Systems with the Computer)¶

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main¶

(Sommersemester 2024)¶

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske¶

Frankfurt am Main 13.01.2024¶

Erster Vorlesungsteil:¶

Eine kleine Einführung in Jupyter Notebooks¶

Einführung¶

Die Programmiersprache Python ist wie C++ eine höhere Objekt-orientierte Programmiersprache. Jedoch finden viele erfahrene Programmierer die Programmiersprache Python weniger strukturiert als C++ und auch, hinsichtlich der Performance der erstellten Programme ist die Programmiersprache C++ ein wenig besser, und große, aufwendige Simulationsprogramme werden oft nicht mittels Python programmiert. Im Bereich der Visualisierung von Daten und im Bereich der Illustration von mathematisch/physikalischen Gleichungen hat Python jedoch sicherlich Vorteile gegenüber C++ und bietet mittels des Python-Moduls "matplotlib" (siehe Matplotlib: Visualization with Python) einen einfachen Weg Bilder und Animationen des simulierten Systems zu erzeugen.

Im Folgenden wird Ihnen keine strukturierte "Einführung in die Programmiersprache Python" gegeben, sondern es werden nur ganz spezielle Themen der Programmiersprache Python vorstellt, die wir in den weiteren Python Jupyter Notebooks verwenden werden. Strukturierte Einführungen in die Programmiersprache Python finden Sie z.B. unter den folgenden Links:

Literatur zu Python¶

  • Python-Onlinekurs auf Deutsch
  • Python 3 documentation
  • Hans Petter Langtangen: A Primer on Scientific Programming with Python
  • David M. Beazley: Python - Essential Reference
  • B. Slatkin: Effective Python

In Python Jupyter Notebooks kann man die Programmiersprache Python in einer anwendungsfreundlichen Umgebung nutzen und die berechneten Ergebnisse auch gleich visualisieren. Sie sind somit in ihrem Erscheinungsbild den kommerziellen Software-Produkten Mathematica oder Maple sehr ähnlich. Die in dieser Vorlesung bereitgestellten Jupyter Notebooks benutzen als Kernel Python 3.

Grundlagen¶

Möchte man eine Gleichung in einem C++ Programm implementieren, muss man zunächst die in ihr enthaltenen Variablen im C++ Programm deklarieren (näheres siehe Datentypen und Variablen in C++). Dies ist in der Programmiersprache Python nicht nötig und beim Initialisierungsprozess einer Variable wird, ähnlich dem C++ auto-Befehl, automatisch der Typ festgelegt.

In [1]:
zeichen = "Hallo"
zahl_1 = 3
zahl_2 = 2.45688586

Den Datentyp einer definierten Variable erhält man mit:

In [2]:
type(zeichen)
Out[2]:
str
In [3]:
type(zahl_1)
Out[3]:
int
In [4]:
type(zahl_2)
Out[4]:
float

Der Wert einer Variable ist unter ihrem Variablennamen gespeichert:

In [5]:
zahl_2
Out[5]:
2.45688586

Eine Zahlenliste in Python wird z.B. wie folgt definiert:

In [6]:
meine_liste = [0,1,4,9,16]
In [7]:
type(meine_liste)
Out[7]:
list

Die Werte der definierten Liste sind unter ihrem Variablennamen gespeichert.

In [8]:
meine_liste
Out[8]:
[0, 1, 4, 9, 16]

Möchte man z.B. den dritten Eintrag der Liste ausgeben lassen, schreibt man

In [9]:
meine_liste[2]
Out[9]:
4

''meine_liste$[2]$'', da der Index der Liste beginnend von Null gezählt wird. Einen weiteren Wert fügt man der Liste mit dem Befehl ''meine_liste.append()'' hinzu.

In [10]:
meine_liste.append(34)
In [11]:
meine_liste
Out[11]:
[0, 1, 4, 9, 16, 34]

Möchte man sich einen speziellen Bereich der Listeneinträge ausgeben lassen (z.B. von 1 bis 4) schreibt man

In [12]:
meine_liste[1:4]
Out[12]:
[1, 4, 9]

Neben Variablen und Listen lassen sich in Python auch Funktionen definieren. Die Definition von Funktionen in Python ist ähnlich der in C++ (näheres siehe Funktionen in C++), jedoch, von der Syntax her, ein wenig einfacher.

Betrachten wir uns z.B. die Funktion $g(x) = x^2 + 2 x + 3$. Eine Definition der Funktion in Python würde wie folgt lauten

In [13]:
def g(x):
    wert = x**2 + 2*x + 3
    return wert

und der Funktionswert an der Stelle $x=2.1$ berechnet man dann mittels $g(2.1)$.

In [14]:
g(2.1)
Out[14]:
11.61

Das Modul NumPy¶

In Python existieren diverse Erweiterungsmodule für spezielle Programierungsanwendungen. Falls Sie Module das erste Mal verwenden, müssen Sie diese zunächst in Ihrer Python 3 Umgebung installieren (z.B. in einem Linux Terminal mit "pip3 install numpy"). Für numerische Berechnungen ist hierbei das Erweiterungsmodul ''NumPy'' sehr hilfreich (näheres siehe z.B. NumPy Tutorial). Bevor wir NumPy benutzen können, müssen wir es importieren. Es ist hierbei üblich, das Modul beim importieren in ''np'' umzubenennen:

In [15]:
import numpy as np

Möchte man nun die im Modul vordefinierten Größen oder Funktionen verwenden, so muss man ''np.'' vor dem jeweiligen Befehl schreiben. Zum Beispiel erhält man den Zahlenwert der Zahl $\pi$ wie folgt:

In [16]:
np.pi
Out[16]:
3.141592653589793

Das Modul NumPy stellt eine nützliche Speichermöglichkeit von nummerischen mehrdimensionalen Arrays (Vektoren, Matrizen,...) bereit, die gerade bei großen Datenmengen Python-Listen vorzuziehen ist, da diese ''NumPy-Arrays'' auf Geschwindigkeit optimiert sind. Mittels der Funktion ''np.linspace(-2, 3, 50)'' kann man z.B. das x-Werteintervall $x \in [-2 , 3]$ in 50 Werte unterteilen und als eindimensionales NumPy-Array darstellen.

In [17]:
x_werte = np.linspace(-2, 3, 50)
In [18]:
x_werte
Out[18]:
array([-2.        , -1.89795918, -1.79591837, -1.69387755, -1.59183673,
       -1.48979592, -1.3877551 , -1.28571429, -1.18367347, -1.08163265,
       -0.97959184, -0.87755102, -0.7755102 , -0.67346939, -0.57142857,
       -0.46938776, -0.36734694, -0.26530612, -0.16326531, -0.06122449,
        0.04081633,  0.14285714,  0.24489796,  0.34693878,  0.44897959,
        0.55102041,  0.65306122,  0.75510204,  0.85714286,  0.95918367,
        1.06122449,  1.16326531,  1.26530612,  1.36734694,  1.46938776,
        1.57142857,  1.67346939,  1.7755102 ,  1.87755102,  1.97959184,
        2.08163265,  2.18367347,  2.28571429,  2.3877551 ,  2.48979592,
        2.59183673,  2.69387755,  2.79591837,  2.89795918,  3.        ])

'NumPy-Arrays'' kann man, falls benötigt, auch wieder in Python-Listen umwandeln: list(x_werte)

Die Funktionswerte der oben bereits definierten Funktion $g(x) = x^2 + 2 x + 3$ in diesem x-Werteintervall, lassen sich nun mit folgendem Befehl berechnen:

In [19]:
g(x_werte)
Out[19]:
array([ 3.        ,  2.8063307 ,  2.63348605,  2.48146606,  2.35027072,
        2.23990004,  2.15035402,  2.08163265,  2.03373594,  2.00666389,
        2.00041649,  2.01499375,  2.05039567,  2.10662224,  2.18367347,
        2.28154935,  2.4002499 ,  2.53977509,  2.70012495,  2.88129946,
        3.08329863,  3.30612245,  3.54977093,  3.81424406,  4.09954186,
        4.40566431,  4.73261141,  5.08038317,  5.44897959,  5.83840067,
        6.2486464 ,  6.67971678,  7.13161183,  7.60433153,  8.09787589,
        8.6122449 ,  9.14743857,  9.70345689, 10.28029988, 10.87796751,
       11.49645981, 12.13577676, 12.79591837, 13.47688463, 14.17867555,
       14.90129113, 15.64473136, 16.40899625, 17.1940858 , 18.        ])

Das Modul Matplotlib¶

Mittels des Python-Moduls Matplotlib (siehe Matplotlib: Visualization with Python) können Grafiken und Animationen erzeugt werden.

In [20]:
import matplotlib.pyplot as plt

Das x-y Diagramm der Funktion $g(x) = x^2 + 2 x + 3$ im Intervall $x \in [-2 , 3]$ lässt sich wie folgt darstellen:

In [21]:
plt.plot(x_werte,g(x_werte))
plt.show()
No description has been provided for this image

Die Visualisierung von Funktionen, die von mehreren Variablen abhängen, kann man auf unterschiedlichen Weisen in Jupyter Notebooks realisieren. Betrachten wir z.B. die Funktion $f(x,y)={\rm sin}(x)+{\rm cos}(y)$

In [22]:
def f(x,y):
    wert = np.sin(x) + np.cos(y)
    return wert

Die Funktion stellt eine Fläche im Raum dar und wir möchten diese nun mittels Matplotlib visualisieren. Hierzu erzeugen wir zunächst einen zweidimensionalen x-y Werte-Raum ($x \in [0 , 12]$ und $y \in [0 , 12]$) mit z.B. $300 \cdot 300$ Punkten.

In [23]:
X = np.linspace(0, 12, 300)
Y = np.linspace(0, 12, 300)
X, Y = np.meshgrid(X, Y)

Für das 3-dimensionale Achsensystem und für die farbliche Höhenvisualisierung der Fläche importieren wir noch zwei weitere Python-Pakete.

In [24]:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

Die Farbe der Fläche soll sich in Abhängigkeit von ihrer Höhe unterscheiden, wobei wir das Farbschema ''Colormap: cm.coolwarm'' wählen, bei dem negative Funktionswerte blau und positive rot erscheinen (näheres siehe Choosing Colormaps). Eine Visualisierung der Fläche und Speicherung des Bildes erhält man z.B. mittels der folgenden Befehle:

In [25]:
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
surf = ax.plot_surface(X, Y, f(X,Y), cmap=cm.coolwarm, linewidth=0, alpha=1)
ax.view_init(azim=35, elev=55)
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("f")
plt.savefig("./Bild.png",format="png")
No description has been provided for this image

Man kann in Matplotlib auch animierte Grafiken erzeugen und diese speichern und in Jupyter Notebooks darstellen (näheres siehe matplotlib.animation und Beispiele von Animationen). Hierzu importieren wir uns die folgenden Pakete.

In [26]:
import matplotlib.animation as animation
from IPython.display import HTML

In der folgenden Animation wollen wir die Fläche der Funktion $h(x,y,t)={\rm sin}(\frac{x}{\frac{t}{25}+1})+{\rm cos}(\frac{y}{\frac{t}{25}+1})$ in einem xyz-Diagramm dargestellen und den Parameter $t$ im Bereich $t \in [0,30]$ variieren. Hierzu definieren wir uns zunächst die Funktion.

In [27]:
def h(x,y,t):
    wert = np.sin(x/(t/25+1)) + np.cos(y/(t/25+1))
    return wert

Im Folgenden wird die Animation erzeugt, als animated GIF gespeichert (Animation.gif) und im Jupyter Notebook dargestellt.

In [28]:
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.view_init(azim=35, elev=55)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("h")
fig.colorbar(surf, shrink=0.5, aspect=10, ax=ax)

def animate(t):
    ax.cla()
    ax.plot_surface(X, Y, h(X, Y,t), cmap=cm.coolwarm, linewidth=0, alpha=1)
    return fig,

ani = animation.FuncAnimation(fig,animate,frames=30,interval=50)

ani.save("./Animation.gif")
plt.close(ani._fig)
HTML(ani.to_html5_video())
Out[28]:
Your browser does not support the video tag.

Die Python Bibliothek Plotly¶

Die Visualisierung kann man auch mittels der "Plotly Python Open Source Graphing Library" erstellen. In dieser Bibliothek ist es möglich, interaktive Abbildungen in Python darzustellen ( siehe Die Python Bibliothek Plotly). Im Folgenden sehen Sie die Visualisierung der Funktion $f(x,y)={\rm sin}(x)+{\rm cos}(y)$ mittels Plotly.

In [29]:
import plotly.graph_objects as go

fig = go.Figure(go.Surface(x=X, y=Y, z=f(X,Y), colorscale='jet', showscale=True))

fig.update_layout(autosize=True,
                  width=700, height=700,
                  margin=dict(l=15, r=50, b=65, t=90))
fig.update_layout(scene = dict(
                    xaxis_title='x',
                    yaxis_title='y',
                    zaxis_title='f'))

Das Modul SymPy¶

Das Modul ''SymPy'' stellt eine umfangreiche Bibliothek für symbolische Mathematik bereit. Dieses Computer-Algebra-System macht symbolische Berechnungen in Python einfach möglich (näheres siehe Das Modul SymPy). Es werden im Folgenden einige Beispiele von symbolischen Berechnungen in Python behandelt.

Das Matrix Berechnungen¶

In [30]:
from sympy import *
init_printing()

Wir definieren uns eine (2x3)-Matrix $\hat{\bf {A}}$ und eine (3x3)-Matrix $\hat{\bf {B}}$:

In [31]:
MA = Matrix([[-1,3,9],[2,-5,6]])
MA
Out[31]:
$\displaystyle \left[\begin{matrix}-1 & 3 & 9\\2 & -5 & 6\end{matrix}\right]$
In [32]:
MB = Matrix([[-7,-1,5],[-9,-3,4],[1,8,-2]])
MB
Out[32]:
$\displaystyle \left[\begin{matrix}-7 & -1 & 5\\-9 & -3 & 4\\1 & 8 & -2\end{matrix}\right]$

Das Produkt der Matrizen $\hat{\bf {A}}$ und $\hat{\bf {B}}$ berechnet sich zu

In [33]:
MA*MB
Out[33]:
$\displaystyle \left[\begin{matrix}-11 & 64 & -11\\37 & 61 & -22\end{matrix}\right]$

, die transponierte Matrix $(\hat{\bf {A}})^T$

In [34]:
transpose(MA)
Out[34]:
$\displaystyle \left[\begin{matrix}-1 & 2\\3 & -5\\9 & 6\end{matrix}\right]$

und die Determinante der Matrix $\hat{\bf {B}}$ berechnet sich zu:

In [35]:
MB.det()
Out[35]:
$\displaystyle -149$

Weitere Matrix Operationen finden Sie unter http://docs.sympy.org/latest/tutorial/matrices.html

Analytische Lösungen von Gleichungen¶

Betrachten wir z.B. die folgende Gleichung $f(x) = x^2-3x+9=0$. Zunächst definieren wir die Variable $x$ als ein Sympy-Symbol und dann die Gleichung:

In [36]:
x = symbols('x')
f = x**4 - 5*x**2 + 2
Gl = Eq( f , 0)
Gl
Out[36]:
$\displaystyle x^{4} - 5 x^{2} + 2 = 0$

Die analytischen Lösungen der Gleichung lauten:

In [37]:
Loes = solve(Gl)
Loes
Out[37]:
$\displaystyle \left[ - \sqrt{\frac{5}{2} - \frac{\sqrt{17}}{2}}, \ \sqrt{\frac{5}{2} - \frac{\sqrt{17}}{2}}, \ - \sqrt{\frac{\sqrt{17}}{2} + \frac{5}{2}}, \ \sqrt{\frac{\sqrt{17}}{2} + \frac{5}{2}}\right]$

Die numerischen Werte erhält man wie folgt:

In [38]:
Loes[0].evalf(),Loes[1].evalf(),Loes[2].evalf(),Loes[3].evalf()
Out[38]:
$\displaystyle \left( -0.662153446861956, \ 0.662153446861956, \ -2.13577920506986, \ 2.13577920506986\right)$

Man kann dies sehen, indem man die linke Seite der Gleichung, die Funktion $f(x)$, mittels Sympy als ein Plot darstellt.

In [39]:
plot(f, (x, -3, 3));
No description has been provided for this image

Differentation und Integration¶

Man kann in Sympy auch Funktionen differenzieren und integrieren. Man bestimmt die erste Ableitung $\frac{d f(x)}{dx}$ einer Funktion mittels ''diff(f,x)'' bzw. ''f.diff(x)''

In [40]:
diff(f,x)
Out[40]:
$\displaystyle 4 x^{3} - 10 x$

und die zweite Ableitung $\frac{d^2 f(x)}{dx^2}$

In [41]:
diff(f,x,x)
Out[41]:
$\displaystyle 2 \left(6 x^{2} - 5\right)$

In ähnlicher Weise bestimmt man das unbestimmte Integral $\int f(x) \,dx$ einer Funktion.

In [42]:
integrate(f,x)
Out[42]:
$\displaystyle \frac{x^{5}}{5} - \frac{5 x^{3}}{3} + 2 x$

Den Wert des bestimmten Integrals erhält man, indem man die Integrationsgrenzen festlegt, z.B. $\int_0^3 f(x) \,dx$:

In [43]:
integrate(f,(x, 0, 3))
Out[43]:
$\displaystyle \frac{48}{5}$

Analytische Lösungen von Differenzialgleichungen¶

Zusätzlich ist es in Sympy auch möglich Differenzialgleichungen (DGLs) zu lösen (falls eine analytische Lösung existiert). Betrachten wir z.B. die folgende DGL zweiter Ordnung: $\frac{d^2y(x)}{dx^2}=-3y(x)$. Zunächst definieren wir die gesuchte Funktion $y(x)$ als ein Sympy-Funktion und dann die zugrundeliegende Differenzialgleichung:

In [44]:
y = Function('y')(x)
DGL = Eq(y.diff(x,x), -3*y)
DGL
Out[44]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = - 3 y{\left(x \right)}$

Die allgemeine Lösung lautet:

In [45]:
dsolve(DGL)
Out[45]:
$\displaystyle y{\left(x \right)} = C_{1} \sin{\left(\sqrt{3} x \right)} + C_{2} \cos{\left(\sqrt{3} x \right)}$

Die allgemeine Lösung hängt von den zwei Parametern $C_1$ und $C_2$ ab, die durch eine Wahl der Anfangs/Randbedingungen festgelegt werden können. Wir wählen z.B.

In [46]:
Eq(y.subs(x,0),3)
Out[46]:
$\displaystyle y{\left(0 \right)} = 3$
In [47]:
Eq(y.diff(x).subs(x, 0),2)
Out[47]:
$\displaystyle \left. \frac{d}{d x} y{\left(x \right)} \right|_{\substack{ x=0 }} = 2$

und erhalten als spezielle Lösung der DGL:

In [48]:
Loes = dsolve(DGL,ics={y.subs(x,0):3,y.diff(x).subs(x, 0): 2})
Loes
Out[48]:
$\displaystyle y{\left(x \right)} = \frac{2 \sqrt{3} \sin{\left(\sqrt{3} x \right)}}{3} + 3 \cos{\left(\sqrt{3} x \right)}$

Die Lösung können wir dann in Sympy wie folgt visualisieren:

In [49]:
plot(Loes.rhs, (x, 0, 10));
No description has been provided for this image

Alternativ können wir die Lösung auch mit Matplotlib darstellen, wobei man hierbei zunächst den analytischen Ausdruck der Lösungsfunktion in eine numerische Funktion umwandelt (''lambdify(x, Loes.rhs)'').

In [50]:
x_werte = np.linspace(0, 10, 100)
Loes_l = lambdify(x, Loes.rhs)
plt.ylabel('y(x)')
plt.xlabel('x')
plt.plot(x_werte,Loes_l(x_werte));
No description has been provided for this image

Weitere Beispiele des Moduls Sympy finden Sie in der Vorlesung Einführung in die Programmierung für Studierende der Physik (Sommersemester 2022) und in dem Jupyter Notebook Das Computeralgebrasystem: Python Jupyter + SymPy