Solving Advanced Problems with Python

Idee & Ziele:

• Verschiedene anspruchsvolle Probleme aus dem mathematisch-naturwissenschaftlichen Bereich mit Python lösen.
• Mit Jupyter-Notebooks arbeiten

Zeit: 1 Woche

1. Studiere das kurze Tutorial zu Jupyter Notebooks und installiere beide Varianten, also Jupyter im Browser und Jupyter in VSCode.
2. Lade das PDF mit dem Jupyter-Cheat Sheet herunter und verschaffe dir eine grobe Übersicht.
3. Lade das folgende Jupyter-Notebook herunter (zuerst entzippen) und arbeite es durch: 01_intro_widgets.ipynb.zip

Zeit: 1 Woche

Erstelle ein Jupyter-NB mit Namen 02_einfache_dgl.ipynb im entsprechenden Repo.

Bei einer Gleichung (hier: quadratische Gleichung) wie z.B. $$x^2 = 20 - x$$ sucht man alle Zahlen, die man für die Variable $x$ einsetzen kann, so dann die Gleichung erfüllt ist. Die Gleichung in diesem Beispiel hat zwei Lösungen: $x=-5$ und $x=4$.

Eine Differentialgleichungen (DGL) ist vom Prinzip her ähnlich. Allerings sucht man als Lösung(en) nicht Zahlen, sondern Funktionen, die eine gewisse Bedingung erfüllen. Ein Beispiel wäre:

$$y'(x) = \frac{2 \cdot y(x)}{x}$$

Hier suchen wir alle Funktionen $y(x)$, deren Ableitung $y'(x)$ gleich ist wie zwei mal die ursprüngliche Funktion dividiert durch die Variable $x$. Die Lösung von dieser Differentialgleichung ist die Funktion

$$y(x) = k \cdot x^2$$

Dabei ist $k$ eine beliebige Zahl. Da wir für $k$ unendlich viele Möglichkeiten haben, gibt es also unendlich viele Lösungen.

Zum Schluss wollen wir uns noch vergewissern, dass die gefundene Funktion tatsächlich die DGL löst:

$$y'(x) = 2 k x = \frac{2kx^2}{x} = \frac{2 \cdot y(x)}{x}$$

Differentialgleichungen von Hand zu lösen kann sehr anspruchsvoll und oft sogar unmöglich sein. Wie dem auch sei, dies ist nicht das, was wir hier machen wollen - dies wird im PHAM behandelt. Hier geht es um das numerische Lösen von Differentialgleichungen, wir überlassen also dem Computer das Lösen.

Wir betrachten die Differentialgleichung

$$y'(x) = - k y(x)$$

wobei $k$ irgend eine konstante Zahl ist. Wir suchen also eine Funktion $y(x)$, die, wenn man sie ableitet, minus sich selbst mal eine fixe Konstante ergibt.

Wir wollen nun diese DGL mithilfe von Python lösen. Dazu verwenden wir den ODEINT Solver. ODE steht für Ordinary Differential Equation (Gewöhnliche Differentialgleichungen) und INT für Integrate. ODEs sind Differentialgleichungen, die von nur einer Variablen, hier $x$, abhängen. Im Gegensatz dazu hängen Partielle Differentialgleichung von mehreren Variablen ab.

Beachte: Es gibt unendlich viele Funktionen $y(x)$, die die DGL erfüllen, alle haben aber die gleiche Form. Um eine eindeutige Lösung zu erhalten, muss man deshalb eine Anfangsbedingung festlegen. Zum Beispiel kann man sagen, dass für $x=0$ der Funktionswert $y=5$ sein soll. Diese Anfangsbedingung muss dem Solver als Argument übergeben werden. Die gesamte Information, die man benötigt, um eine eindeutige Lösung zu finden, sieht dann wie folgt aus:

$$y'(x) = - k y(x)$$ $$y(0) = 5$$

Wichtig ist, dass in der DGL die Ableitung alleine auf der linken Seite steht.

Links:

• Manual ODEINT: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

Lösung mit odeint:

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
 
def dgl(y,x):
    k = 0.4
    y_prime = - k * y
    return y_prime
 
X = np.linspace(0,20,100)
 
y0 = 15
 
Y = odeint(dgl,y0,X)
 
plt.plot(X,Y)

• Wissen, was DGL sind …
• und wozu sie nützlich sind
• Wissen, warum DGL Anfangsbedingungen brauchen
• DGL mit odeint lösen können
• Mit matplotlib ansprechende Figuren erstellen können

1. Gehe das Beispiel oben Schritt für Schritt durch, stelle sicher, dass du alles verstehst. Studiere dazu z.B. das Manual von ODEINT.
2. Um was für eine Funktion $y(x)$ handelt es sich hier?
3. Spiele mit dem Wert für den Parameter $k$. Welchen Einfluss hat dieser auf die Funktion?
4. Was kann im real life durch eine solche Funktion beschrieben werden? Mache ein Beispiel

1. Verschönere die Figur, unter anderem sollst du folgendes machen:
   1. Achsen beschriften
   2. Kurve dicker und andere Farbe
   3. Legende hinzufügen
   4. Titel hinzufügen
2. Nun soll die DGL für verschiedene Werte für $k$ gelöst werden. Finde dazu einen Weg, wie dem odeint-Solver ein zusätzliches Argument übergeben werden kann. Visualisiere dann alle Funktionsgraphen in einer Figur und beschrifte die Graphen in der Legende.

Mache eine Kopie des Codes und füge einen Slider ein, mit dem man den Wert von $k$ einstellen kann.

1. Löse die Differentialgleichung $$f'(x) = a$$ mit odeint, wobei $a$ ein beliebiger Parameter ist.
2. Stelle die Funktion $f(x)$ graphisch dar.
3. Um was für eine Funktion handelt es sich bei $f(x)$?

1. Löse die DGL vom Theorieteil am Anfang des Kapitels: $$f'(x) = \frac{2 \cdot f(x)}{x}$$
2. Stelle die Lösung graphisch dar und vergleiche sie mit der exakten Lösung $f(x) = k x^2$.

Finde heraus, wie du die Darstellung deiner Figur verändern kannst, z.B.:

1. Grösse der Figur
2. Farbe und Liniendicke von Graph
3. Hintergrundfarbe
4. Ticks auf Achsen
5. Beschriftungen: Titel, Achsen, Graph, Legende
6. …

Zeit: 1 Woche

In der Physik beschreibt man die Bewegungen von Objekten (Beispiele: Stein fällt nach unten, Masse an Feder die oszilliert, …) durch sogenannte Bewegungsgleichungen, welche Differentialgleichungen (DGL) sind.

Dabei ist es wichtig zu wissen, wie Position, Geschwindigkeit und Beschleunigung zusammenhängen:

Sei $x(t)$ die Position eines Objekts als Funktion der Zeit. Die Geschwindigkeit des Objekts zu jedem Zeitpunkt ist dann gegeben durch die Ableitung der Position: $$v(t) = x'(t)$$ Die Geschwindigkeit gibt also an, wie fest sich die Position mit der Zeit verändert. Ganz analog ist die Beschleunigung die Änderung der Geschwindigkeit und ist deshalb gegeben durch $$a(t) = v'(t) = x''(t)$$ Die Beschleunigung ist also die zweite Ableitung der Position.

Aus dem Physikunterricht kennst du sicher das Hookesche Gesetz: $$ F = - k x $$ Es beschreibt die Bewegung eines Federpendels ohne Reibung. Das Hookesche Gesetz gibt uns die Federkraft $F$ eines Federpendels an. Es besagt, dass die Federkraft gegeben ist durch minus die Federkonstante $k$ mal die Auslenkung der Feder aus ihrer Gleichgewichtsposition. Das Minus sagt, dass die Kraft immer gegen der Auslenkung wirkt: Ist die Feder gedehnt, so will die Kraft die Feder zusammendrücken. Ist sie aber zusammengedrückt, will die Kraft sie wieder auseinander ziehen.

Das Hookesche Gesetz können wir nun als Differentialgleichung schreiben. Die Beschleunigung der Feder $a$ und die Federkraft $F$ hängen über das zweite Newtonsche Gesetz zusammen: $F = m a$, wobei $m$ die Masse des Objekts ist, welches an der Feder befestigt ist. Damit können wir das Hookesche Gesetz wie folgt umschreiben:

$$F = - k x$$

$$m a = - k x$$

$$m x''(t) = - k x(t)$$

$$x''(t) = - \frac{k}{m} x(t)$$

Wir haben nun also wieder eine Differenzialgleichung. Im Gegensatz zu den DGL von der letzten Übung enthält diese aber eine zweite Ableitung. Zur Erinnerung: Die DGL von der letzten Übung, welche erste Ableitungen beinhaltet haben, hatten immer unendlich viele Lösungen. Wir mussten deshalb eine Randbedingung festlegen.

Löst man nun die neue DGL (Hook'sches Gesetz), so erhält man die Position der Feder als Funktion der Zeit $x(t)$. Hier gibt es wieder unendlich viele Möglichkeiten. Da die neue DGL eine zweite Ableitung beinhaltet, muss man nun zwei Randbedingungen angeben - eine für $x(t)$ und eine für deren Ableitung $x'(t)$, z.B.:

$$x(0) = x_0$$

$$x'(0) = 0$$

Man muss also einen Zeitpunkt (typischerweise den Anfang) wählen und sagen wo und wie schnell das Federpendel dann ist. Wir legen also fest, dass die Feder zum anfänglichen Zeitpunkt $t=0$ an der Position $x_0$ ist und die Geschwindigkeit dann $=0$ ist. Dies macht Sinn, wenn man von Hand eine Feder spannt und zum Zeitpunkt $t=0$ los lässt.

Diese Differentialgleichung kann man analytisch, also exakt, lösen. Mit den obigen Randbedingungen erhält man

$$x(t) = x_0 \cos\left(\sqrt{\frac{k}{m}} \cdot t\right)$$

Ein System, welches so schwingt, nennt man einen harmonischen Oszillator und seine Bewegung eine harmonische Bewegung. Diese Bewegung ist ungedämpft, das Federpendel schwingt also unendlich lange.

Die DGL für den Harmonischen Oszillator enthält eine zweite Ableitung. Wie können wir solche mit odeint lösen? Der Trick ist folgender: Wandle die eine DGL mit einer zweiten Ableitung um in zwei DGL mit nur ersten Ableitungen:

Anstelle von

$$x''(t) = - \frac{k}{m} x(t)$$

betrachten wir

$$x'(t) = v(t)$$

$$v'(t) = - \frac{k}{m} x(t)$$

Nun können wir diese beiden DGL zusammen mit odeint lösen. Die Funktion, die odeint übergeben wird, sieht dann wie folgt aus:

def harmonic(x_v,t,k,m):
    x = x_v[0]
    v = x_v[1] 
    return [v, -k/m * x]

Beachte, dass x_v eine Liste ist. Das erste Element beinhaltet die Position $x$, das zweite die Geschwindigkeit $v$ (also die Ableitung von $x$).

Der harmonische Oszillator ist gut und schön, in Realität tritt aber meist Reibung auf, welche die Bewegung dämpft. Die Amplitude (Auslenkung) eines Federpendels nimmt aufgrund von Reibungskräften wie Luftwiderstand kontinuierlich ab.

Wir können in unserer DGL nun Reibung ganz einfach einbauen, in dem wir einen entsprechende Kraftterm in unsere DGL einbauen. Zum Beispiel hängt der Luftwiderstand von der Geschwindigkeit ab: Je höher die Geschwindigkeit, desto grösser der Luftwiderstand. Wir führen deshalb einen Reibungsterm ein, der linear ist in der Geschwindigkeit: $$F = - k x - c v$$ Dabei ist $c$ eine Konstante, die uns angibt, wie stark der Luftwiderstand ist.

Erstelle ein Jupyter-NB mit Namen 03_bewegungsgleichungen.ipynb im entsprechenden Repo.

1. Beweise, dass $x(t) = x_0 \cos\left(\sqrt{\frac{k}{m}} \cdot t\right)$ eine Lösung des harmonischen Oszillators ist: Setze einfach in die DGL ein und zeige, dass links und rechts das Gleiche herauskommt.
2. Simuliere nun einen harmonischen Oszillator $x''(t) = - \frac{k}{m} x(t)$ wie ein Federpendel ohne Reibung, in dem du die DGL mit odeint löst.
3. Stelle die Lösungsfunktion $x(t)$ dann mit matplotlib schön graphisch dar.
4. Vergleiche nun deine numerische Lösung mit der analytischen Lösung $x(t) = x_0 \cos\left(\sqrt{\frac{k}{m}} \cdot t\right)$ und stelle sicher, dass diese übereinstimmen.

1. Simuliere nun den gedämpften harmonischen Oszillator mit einem Reibungsterm, der linear von der Geschwindigkeit abhängt.
2. Füge nun Slider (widgets) hinzu, mit denen du die Werte für die Konstanten $k,m,c$, die Anfangsposition $x_0$, sowie die Zeit-Range einstellen kannst. Sobald ein Slider geändert wurde, soll der Graph entsprechend angepasst werden. Siehe Tipps unten.

Tipps für Widgets:

• Verwende widgets.interactive(update_plot,param1,param2,...)
• definiere Funktion update_plot(param1,param2,...), die Lösung der DGL und Darstellung der Lösung enthält
• param1 usw. sind die Parameter, die mit Widgets (z.B. Slider) eingestellt werden

Zeit: 2 Wochen

Gegeben sind zwei Massen $m_1$ und $m_2$ mit …

• Anfangsposition $\vec{r}_1(t=0)$ und $\vec{r}_2(t=0)$
• Anfangsgeschwindigkeit $\vec{v}_1(t=0)$ und $\vec{v}_2(t=0)$

Ziel: Bestimme Position der beiden Massen zu jedem Zeitpunkt, also $\vec{r}_1(t) = ?$ und $\vec{r}_2(t) = ?$

Die beiden Massen ziehen sich natürlich gegenseitig mit der Gravitationskraft an: $$F_G = G \frac{m_1 \cdot m_2}{d^2}$$ wobei $d$ der Abstand zwischen den beiden Massen ist.

Da eine Kraft sowohl eine Stärke als auch eine Richtung hat, wird sie als Vektor dargestellt. Die Formel oben, die du bereits aus dem GF Physik kennst, gibt allerdings nur die Stärke an. Wir brauchen aber die Vektor-Version: Die Gravitationskraft, die wegen $m_2$ auf $m_1$ wirkt, ist: $$\vec{F}_{12} = G \frac{m_1 \cdot m_2}{|\vec{r}_{12}|^3} \vec{r}_{12}$$ wobei $\vec{r}_{12} = \vec{r}_2 - \vec{r}_1$ der Verschiebungsvektor ist, der von $m_1$ zu $m_2$ zeigt. Damit ergibt sich: $$\vec{F}_{12} = G \frac{m_1 \cdot m_2}{|\vec{r}_2 - \vec{r}_1|^3} (\vec{r}_2 - \vec{r}_1)$$

Wahrscheinlich fragst du dich, woher das „hoch 3“ kommt! Dies macht Sinn, wenn man die Stärke dieses Vektors bestimmt: $$|\vec{F}_{12}| = G \frac{m_1 \cdot m_2}{|\vec{r}_{12}|^3} |\vec{r}_{12}| = G \frac{m_1 \cdot m_2}{|\vec{r}_{12}|^2}$$

Die Differenzialgleichungen (DGL), welche die Position der Masse $m_1$ bestimmt ist dann wie immer durch das 2. Newtonsche Gesetz gegeben: $$m_1 \, \ddot{\vec{r}}_1 = \vec{F}_{12} = G \frac{m_1 \cdot m_2}{|\vec{r}_2 - \vec{r}_1|^3} (\vec{r}_2 - \vec{r}_1)$$ und gleichermassen für die zweite Masse: $$m_2 \, \ddot{\vec{r}}_2 = \vec{F}_{21} = G \frac{m_1 \cdot m_2}{|\vec{r}_1 - \vec{r}_2|^3} (\vec{r}_1 - \vec{r}_2)$$

In einem Zweikörpersystem bewegen sich die beiden Massen immer in einer Ebene. Wir können daher das Problem in 2D anstelle 3D betrachten. Die beiden Positionsvektoren haben dann also je eine $x-$ und $y-$ Komponente: $$\vec{r}_1 = \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} \,, \quad \vec{r}_2 = \begin{pmatrix} x_2 \\ y_2 \end{pmatrix}$$ Die DGL sind dann also: $$m_1 \begin{pmatrix} \ddot{x}_1 \\ \ddot{y}_1 \end{pmatrix} = G \frac{m_1 \cdot m_2}{d^3} \begin{pmatrix} x_2 - x_1 \\ y_2 - y_1 \end{pmatrix}$$ $$m_2 \begin{pmatrix} \ddot{x}_2 \\ \ddot{y}_2 \end{pmatrix} = G \frac{m_1 \cdot m_2}{d^3} \begin{pmatrix} x_1 - x_2 \\ y_1 - y_2 \end{pmatrix}$$ $$d = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}$$ Beachte, dass es sich um zwei Vektor-DGL handelt, es sind insgesamt also $4$ DGL: $$\ddot{x}_1 = \frac{G \cdot m_2}{d^3} (x_2 - x_1)$$ $$\ddot{y}_1 = \frac{G \cdot m_2}{d^3} (y_2 - y_1)$$ $$\ddot{x}_2 = \frac{G \cdot m_1}{d^3} (x_1 - x_2)$$ $$\ddot{y}_2 = \frac{G \cdot m_1}{d^3} (y_1 - y_2)$$ Löst man nun diese DGL (z.B. mit odeint), so erhält man die Position der beiden Massen.

Der Schwerpunkt (barycenter) der beiden Massen ist gegeben durch: $$\vec{R} = \frac{m_1\,\vec{r}_1 + m_2\,\vec{r}_2}{m_1 + m_2}$$

Wie oben besprochen geht es im Zweikörperproblem darum, die zu DGL lösen, welche durch das 2. Newtonsche Gesetz gegeben sind:

$$m_1\,\ddot{\vec{r}}_1 = \vec{F_{12}} = \text{Kraft, die Masse 2 auf Masse 1 ausübt}$$ $$m_2\,\ddot{\vec{r}}_2 = \vec{F_{21}} = \text{Kraft, die Masse 1 auf Masse 2 ausübt}$$

Dase 3. Newtonsche Gesetz besagt, dass $$\vec{F_1} = - \vec{F_2}$$ und damit $$\vec{F_1} + \vec{F_2} = 0$$ Es gilt also: $$m_1\,\ddot{\vec{r}}_1 + m_2\,\ddot{\vec{r}}_2 = 0 \quad | \quad :(m_1+m_2)$$ $$\frac{m_1\,\ddot{\vec{r}}_1 + m_2\,\ddot{\vec{r}}_2}{m_1 + m_2} = 0$$ $$\ddot{\vec{R}} = 0$$ Das bedeutet also, dass der Schwerpunkt des 2-Körpersystems nicht beschleunigt ist. Seine Geschwindigkeit $\vec{v} = \dot{\vec{R}}$ ist also konstant.

Wenn man die Lösung des Zweikörperproblems graphisch darstellt, so macht es Sinn, dafür ein Bezugssystem zu wählen, in welchem der sich der Schwerpunkt in Ruhe. Dafür geht man wie folgt vor:

1. Löse die DGL, welche für die beiden Massen durch die Newtonsche Gravitationskraft gegeben sind, ganz normal mit odeint.
2. Berechne dann mithilfe der erhaltenen Lösung für jeden Zeitpunkt den Schwerpunkt des Zweikörpersystems.
3. Subtrahiere den Schwerpunkt von den Positionen der beiden Massen.
4. Plotte die Bahnen der beiden Massen. Du solltest jetzt sehen, dass sie sich auf Ellipsen bewegen.

1. Simuliere die Umlaufbahnen zweier Planeten im Zweikörpersystem. Rechne den Schwerpunkt aus den $x$- und $y$-Koordinaten heraus. Verwende folgende Werte:
   1. $G=1$
   2. Massen: z.B. in Grössenordnung $1-10$.
   3. Startposition: ebenfalls $1-10$.
   4. Anfangsgeschwindigkeit eher $0.05-0.5$
2. Finde dann durch ausprobieren möglichst unterschiedliche Orbits. U.a. sollen die folgenden reproduziert werden:
3. Zweikörpersimulator mit Widgets: Füge einige Widgets ein, mit denen man den Orbit festlegen kann. Sinn würde z.B. folgendes machen:
   1. Massen m1 und m2
   2. Distanz $d$: m1 hat dann Anfangsposition (-d/2,0) und m2 hat dann Anfangsposition (d/2,0).
   3. Anfangsgeschwindigkeiten der beiden, aber immer in $y$-Richtung.
4. (Optional) Animation: Umlaufbahnen sollen animiert dargestellt werden

Zeit: 1 Woche

• Studiere nur Spielregeln
• Implementiere das Game of Life mit Python (normal oder Jupyter)
• Wie du es visualisierst, ist dir überlassen: Konsole, PyGame, PyQt, Jupyter, …
• Programmiere ganz alleine (keine Inspiration von Google, weiter unten auf Wiki, …)

Zeit: 1 Woche

Vorgaben Version v02

Theorie Model vs. View

Empfohlene Klassen

# MODELL
"""
Die ersten beiden Klassen bilden das Modell und haben nichts mit der Darstellung zu tun.
Deshalb dürfen diese KEINE print, PyGame usw. Befehle beinhalten.
"""
class World:
    """
    Create only one object of this class.
    Contains list with all cells that are alive, updated them, ...
    """
    def __init__(self,...):
        self.cells = [] # saves all cells that are ALIVE
        ...
 
    def update(self):
        pass
 
    ....
 
class Cell:
    """
    describes a single cell
    """
    def __init__(self,...):
        self.x = ... # x coordinate
        self.y = ... # y coordinate
        self.is_alive = ... # is currently alive
        self.stays_alive = ... # will still be alive in next round
        self.neighbours = []
 
    ....
 
# VIEW
"""
Die Klasse(n) hier ist für die Darstellung zuständig. 
"""
 
class ConsoleApp: # or PyGameApp
    pass
 
# EXECUTE PROGRAM
"""
Könnte etwa wie folgt aussehen:
"""
 
if __name__ == "__main__":
    world = World("static 1") # erstelle Welt mit vordefiniertem Preset, weitere Möglichkeiten: "glider 1", "random", ...
    console_app = ConsoleApp(world) # erstelle eine Console App (oder PyGame, ...)
    console_app.animate() # führe aus

Zeit: 1 Woche

Ziel: Identischen Model-Code von V2 verwenden für mind. zwei verschiedene Views: Console (aus V2) und GUI (z.B. PyGame/PyQT/Matplotlib)

1. Verbessere falls nötig Code von V2: Model- und View-Code müssen strikt getrennt sein und alles muss strikt objektorientiert programmiert sein.
2. Strukturiere nun Code um. Je ein eigenes File für:
   1. Ganzen Model-Code, z.B. game_of_life_model.py
   2. Code für Console-View, z.B. game_of_life_console.py
   3. Code für GUI-View, z.B. game_of_life_pygame.py
3. Importiere in den beiden Files für die View den Model-Code, z.B. from game_of_life_model import *
4. Programmiere nun die beiden Views aus.
5. Optional: Coole Features einbauen

Erstelle im Eiltempo (2L und ein bisschen zuhause) ein einfaches PyGame-Game, welches an Moorhuhnjagd erinnern soll:

• Statisches Hintergrundbild mit Winterlandschaft
• Es werden zufällig Schneeflocken am oberen Rand erzeugt, die in unterschiedlichen Tempi herunter fallen.
• Es werden zufällig Geschenke am oberen Rand erzeugt, die in unterschiedlichen Tempi herunter fallen.
• Durch Mausklick soll man diese einfangen können: Verschwinden bei Klick, Score erhöht sich um 1.

Weitere optionale Idee:

• Verschiedene Geschenke, gewisse Typen beinhalten Bomben → Abzug oder Leben verlieren bis Game Over
• Santa der mit Schlitten quer durch den Bildschirm reitet
• Power-Ups
• Eigene Ideen, sei kreativ.

Anforderungen Programmieren:

• Muss strikt objektorientiert sein. Mache z.B. Klasse Item sowohl für Geschenke als auch Schneeflocken. Allerdings sollen nur die Geschenke einsammelbar sein. Mache dazu z.B. eine Flag is_collectable, welche True für Geschenke und False für Schneeflocken ist.
• Tipp: Starte mit Template von hier.
• Verwende z.B. diese Bilder oder finde/erzeuge eigene: