Solving Advanced Problems with Python [KSR Wiki]

Diese Seite ist nicht editierbar. Du kannst den Quelltext sehen, jedoch nicht verändern. Kontaktiere den Administrator, wenn du glaubst, dass hier ein Fehler vorliegt.
====== - Solving Advanced Problems with Python ======

**Idee & Ziele:**

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



===== Nützliche Links =====

   * [[pham_physik_sca:python|NumPy, MatPlotLib, Jupyter, Widgets usw. für PHAM]]

===== - Einfache Differentialgleichungen (DGL) =====

Zeit: 1 Woche
==== Einführung ====

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.

==== DGL lösen mit odeint ====

=== Beispiel ===

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:
<code python>
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

def dgl(y,x): # contains right side of 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)
</code>

==== Lernziele ====

   * 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


==== Aufgaben A ====

=== Aufgabe A0 ===

   1. Erstelle auf GitHub ein Repo für diesen Kurs, z.B. `advanced_python_problems`, klone es auf deinen Computer und gib es der LP frei.
   1. Erstelle ein File mit Namen `01_einfache_dgl.ipynb` im entsprechenden Repo.
   1. Richte Jupyter Notebooks in VSCode ein (falls noch nicht gemacht)

=== Aufgabe A1: Beispiel oben ===

== Part I ==

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

== Part II ==

   1. Verschönere die Figur, unter anderem sollst du folgendes machen:
       1. Achsen beschriften
       1. Kurve dicker und andere Farbe
       1. Legende hinzufügen
       1. Titel hinzufügen
   1. 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.

== Part III ==

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


=== Aufgabe A2 ===

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

=== Aufgabe A3 ===

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

=== Aufgabe A4: Matplotlib ===

Gehe zurück zum Plot der letzten Aufgabe und verschönere diese, z.B.:

   1. Grösse der Figur
   1. Farbe und Liniendicke von Graph
   1. Hintergrundfarbe
   1. Ticks auf Achsen
   1. Beschriftungen: Titel, Achsen, Graph, Legende
   1. ...



===== - Bewegungsgleichungen =====

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.

==== Theorie ====

=== Position, Geschwindigkeit & Beschleunigung ===

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.

== Harmonischer Oszillator ==

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.

== DGL mit 2. Ableitungen lösen ==

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:

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

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

== Gedämpfter harmonischer Oszillator ==

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$$
$$mx'' = - k x - c x'$$
$$x'' = - \frac{k}{m} x - \frac{c}{m} x'$$
Dabei ist $c$ eine Konstante, die uns angibt, wie stark der Luftwiderstand ist.

==== Aufgaben B ====

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

=== Aufgabe B1: Harmonischer Oszillator ===

   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.

=== Aufgabe B2: Gedämpfter harmonischer Oszillator ===

   1. Simuliere nun den gedämpften harmonischen Oszillator mit einem Reibungsterm, der linear von der Geschwindigkeit abhängt.
   1. Füge nun Slider (widgets, siehe Link am Anfang) 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.

===== - Chaos: Das Doppelpendel =====

Zeit: 1 Woche

Stell dir vor, du beobachtest ein einfaches Pendel – wie eine Schaukel. Es schwingt regelmässig hin und her, und du kannst ziemlich genau vorhersagen, wo es in einer Sekunde sein wird. Aber was passiert, wenn du ein zweites Pendel an das erste hängst? Dann erhältst du ein Doppelpendels, welches chaotisches Verhalten an den Tag legt.

Ein **Doppelpendel** besteht also aus zwei Pendeln, die miteinander verbunden sind: Das erste Pendel hängt an einem festen Punkt, und das zweite Pendel hängt am Ende des ersten. Wenn du es loslässt, beginnt es sich zu bewegen – aber nicht so ruhig und vorhersehbar wie ein einzelnes Pendel. Stattdessen sieht es oft wild, unregelmässig und fast zufällig aus, es verhält sich *chaotisch*.

In der Physik bedeutet **'chaotisch'** nicht einfach 'durcheinander' oder 'zufällig'. Es bedeutet, dass ein System *extrem empfindlich auf Anfangsbedingungen* reagiert. Wenn du das Doppelpendel also nur ein *winziges bisschen anders startest* – z.B. einen Millimeter höher – kann sich der gesamte Bewegungsverlauf komplett verändern.
Dieses Verhalten nennt man auch den **Schmetterlingseffekt**: Ein kleiner Unterschied am Anfang kann riesige Auswirkungen haben – wie der berühmte Flügelschlag eines Schmetterlings, der einen Sturm auslöst.

Das Doppelpendel zeigt uns, dass selbst einfache physikalische Systeme unvorhersehbar werden können, obwohl sie vollständig durch Naturgesetze beschrieben sind. Es ist ein perfektes Beispiel für **deterministisches Chaos**: Die Regeln sind klar, aber das Ergebnis ist schwer vorherzusagen.

Ziel hier ist, das Doppelpendel zu simulieren: Wir integrieren also seine Bewegungsgleichungen (Differenzialgleichungen) und erstellen damit eine schöne Animation.

{{ :talit:double-pendulum.png?200 |}}


=== Auftrag ===

In Jupiter-Notebook: `doppelpendel.ipynb`

   1. **Integriere die Bewegungsgleichungen** des Doppelpendels mit odeint, genau so, wie du es bereits vom Harmonischen Oszillator kennst. Stelle die beiden Winkel $\theta\_1$ und $\theta\_2$ als Funktionen der Zeit dar. Verzichte hier auf AI.
   1. **Überprüfe**, ob alles korrekt ist, in dem du zu jedem Zeitpunkt die **Gesamtenergie** $E\_\text{tot} = E\_\text{kin} + E\_\text{pot}$ berechnest. Diese muss aufgrund der wohl bekannten Energieerhaltung konstant bleiben.
   1. Erstelle nun eine schöne **Simulation** des Pendels, z.B. mit PyGame. Hier ist AI-Hilfe erlaubt.
   1. Sorge dafür, dass die Animation beliebig lange läuft: Integriere über einen relativ kurzen Abschnitt, integriere danach erneut, wobei du die letzten Werte vom vorherigen Abschnitt als Anfangswerte verwendest.
   1. Optional: Erstelle ein artistisches Video des Doppelpendels mithilfe eines Programms, welches beliebig viele Bilder überlagern kann.


=== Gleichungen und Formeln ===

Die Position der beiden Massen wird vollständig durch die beiden Winkel $\theta\_1$ ('Theta') und $\theta\_2$ bestimmt. Deshalb sucht man nach Bewegungsgleichungen für diese beiden Winkel. Diese beschreiben also, wie sich diese beiden Winkel mit der Zeit verändern. 

**Bewegungsgleichungen:**
$$\ddot{\theta_1} = \frac{ -g\cdot (2\cdot m_1+m_2)\cdot \sin(\theta_1) - m_2\cdot g\cdot \sin(\theta_1-2\cdot \theta_2) - 2\cdot \sin(\theta_1-\theta_2)\cdot m_2\cdot (\dot{\theta_2}^2\cdot l_2 + \dot{\theta_1}^2\cdot l_1\cdot \cos(\theta_1-\theta_2)) }{ l_1\cdot (2\cdot m_1+m_2-m_2\cdot \cos(2\cdot \theta_1-2\cdot \theta_2))}$$
$$\ddot{\theta_2} = \frac{ 2\cdot \sin(\theta_1-\theta_2) \cdot  ( \dot{\theta_1}^2\cdot l_1\cdot (m_1+m_2) + g\cdot (m_1+m_2)\cdot \cos(\theta_1) + \dot{\theta_2}^2\cdot l_2\cdot m_2\cdot \cos(\theta_1-\theta_2) ) }{ l_2\cdot (2\cdot m_1+m_2-m_2\cdot \cos(2\cdot \theta_1-2\cdot \theta_2))}$$

**Koordinaten:**

$$x_1 = l_1 \cdot  \sin(\theta_1)$$
$$y_1 = - l_1 \cdot  \cos(\theta_1)$$
$$x_2 = x_1 + l_2 \cdot  \sin(\theta_2)$$
$$y_2 = y_1 - l_2 \cdot  \cos(\theta_2)$$

**Geschwindigkeit:**
$$v_1 = l_1 \cdot \dot{\theta_1}$$
$$v_2 = \sqrt{l_1^2\cdot \dot{\theta_1}^2 + l_2^2\cdot \dot{\theta_2}^2 + 2\cdot l_1\cdot l_2\cdot \dot{\theta_1} \cdot \dot{\theta_2} \cdot \cos(\theta_1-\theta_2) }$$

**Energie:**
$$E_\text{kin} = \frac12\cdot m_1\cdot v_1^2 + \frac12 \cdot m_2 \cdot v_2^2$$
$$E_\text{pot} = -m_1\cdot g\cdot l_1\cdot \cos(\theta_1) - m_2\cdot g\cdot (l_1\cdot \cos(\theta_1) + l_2\cdot \cos(\theta_2))$$
$$E_\text{tot} = E_\text{kin} + E_\text{pot}$$

**Quellen:**

   * [[https://web.mit.edu/jorloff/www/chaosTalk/double-pendulum/double-pendulum-en.html]]
   * [[https://www.engineered-mind.com/engineering/double-pendulum-1/#implementationinmatlab]]


===== - Chaos II: Der Lorenz-Attraktor =====

Der Lorenz‑Attraktor Der Lorenz‑Attraktor ist ein einfaches deterministisches System von drei gewöhnlichen Differentialgleichungen:
$$\dot{x} = \sigma (y - x)$$
$$\dot{y} = x(\rho - z) - y$$
$$\dot{z} = xy - \beta z$$
Er dient als Prototyp für chaotisches Verhalten. Edward N. Lorenz (MIT, 1963) entwickelte das Modell als stark vereinfachtes Modell für atmosphärische Konvektion. Trotz seiner einfachen Form zeigt das System komplexe, nichtperiodische Bahnen mit empfindlicher Abhängigkeit von den Anfangsbedingungen („Schmetterlingseffekt“).

Typische Parameter (klassisch) $\sigma = 10, \rho = 28, \beta = 8/3$ — diese Wahl führt zum bekannten chaotischen Attraktor.

{{ :talit:lorenz_attractor.png?400 |}}

Programmiere den Lorenz-Attraktor und stelle das Bild schön dar. Ob das Bild auf 1x oder fortlaufend angezeigt wird, ist dir überlassen. Erstelle dazu ein Jupiter-Notebook: `lorenz_attraktor.ipynb`

===== - Eigener DGL Solver =====

Wir haben bereits gelernt, wie man mithilfe von odeint Differenzialgleichungen (DGL) numerisch lösen kann. Nun wollen wir **eigene DGL-Solver** programmieren. Dadurch verstehen wir, was odeint im Hintergrund eigentlich (so ungefähr) macht.

Dazu betrachten wir das folgende **Anfangswert-Problem**:

$$y'(t) = f(t,y(t))$$

$$y(t_0) = y_0$$

Es beinhaltet also eine gewöhnliche **Differentialgleichung erster Ordnung** (nur erste Ableitung) sowie eine **Anfangsbedingung.**
Anschliessen werden verschiedene Lösungsverfahren erläutert:

=== Explizites Euler-Verfahren ===

[[https://de.wikipedia.org/wiki/Explizites_Euler-Verfahren#Das_Verfahren]]

Das einfachste (und auch schlechteste) Verfahren ist das **explizite Euler-Verfahren**.

Man wählt eine diskrete Schrittweite $h>0$ und betrachte dann die Zeitpunkte

$$t_k = t_0 + k \cdot h$$

Die zugehörigen Funktionswerte der Lösungsfunktion sind dann gegeben durch

$$y_{k+1} = y_k + h \cdot f(t_k,y_k)$$

**Achtung:** Wie bei allen numerischen Lösungsverfahren ist die Lösung nicht zu 100\% exakt. Die Werte $y_{k+1}$ weichen von der tatsächlichen, mathematisch korrekten Lösung $y(t_k)$ ab. Je kleiner die Schrittweite $h$ gewählt wird, desto genauer aber stimmen sie überein.


=== Verbessertes explizites Euler-Verfahren ===

[[https://de.wikipedia.org/wiki/Explizites_Euler-Verfahren#Verbessertes_explizites_Euler-Verfahren]]

Mit einer einfachen Modifikation des expliziten Euler-Verfahrens erhält man bessere Resultate:

$$y_{k+\frac12} = y_k + \frac{h}{2} f(t_k,y_k)$$

$$y_{k+1} = y_k + h \cdot f\left(t_k+\frac{h}{2},y_{k+\frac12}\right)$$

=== Runge-Kutta ===

[[https://de.wikipedia.org/wiki/Runge-Kutta-Verfahren#Beispiel]]

Wie beim Euler-Verfahren oben, wählt man eine Schrittweite. Allerdings wählt man nun auch eine beliebige Anzahl von gewichteten Zwischenstufen. Das **dreistufige Runge-Kutta Verfahren** sieht dann wie folgt aus: 

$$y_{k+1} = y_k + h \cdot \left(\frac16k_1 + \frac46k_2 + \frac16k_3 \right)$$

Wobei:

$$k_1 = f(t_k,y_k)$$

$$k_2 = f(t_k+\frac h2,y_k + \frac h2 k_1)$$

$$k_1 = f(t_k+h,y_k-h\cdot k_1 + 2h\cdot k_2)$$

==== Aufträge ====

In Jupiter-Notebook: `eigener_dgl_solver.ipynb`

=== Auftrag 1 ===

Implementiere die drei oben vorgestellten **Lösungsverfahren** für **gewöhnliche Differentialgleichungen erster Ordnung**. Wende sie auf Beispiele aus früheren Lektionen an. Achtung: Bewegungsgleichungen (beinhalten 2. Ableitungen) können damit noch nicht gelöst werden.

Vergleiche dann die Resultate dieser drei Verfahren sowie der odeint-Funktion miteinander. Stelle dazu alle 4 Lösungen gemeinsam in einem Graph dar.

=== Auftrag 2 ===

Modifiziere eines der Lösungsverfahren von oben so, dass du Differentialgleichungen 2. Ordnung, also mit 2. Ableitungen, lösen kannst. Löse damit z.B. den harmonischen Oszillator.

===== - Gravitationspotenzial =====

**Gravitationspotenzial einer homogenen Kugel** (Masse $m_k$, Radius $R_k$ an Position $(x_k,y_k)$):
$$\Phi_k(x,y) = - \frac{G m_k}{r} \quad \text{falls } r > R_k$$
$$\Phi_k(x,y) = - \frac{G m_k}{2\cdot R_k^3}(3\cdot R_k^2 - r^2) \quad \text{falls } r \leq R_k$$
wobei $r = \sqrt{(x-x_k)^2 + (y-y_k)^2}$.

**Gravitationspotential von $n$ homogenen Kugeln**:
$$\Phi(x,y) = \sum_{k=1}^{n} \Phi_k(x,y)$$

Da das Gravitationspotenzial an vielen Punkten berechnet werden soll, macht es Sinn **[[talit:python_advanced#numba|Numba]]** zu verwenden, um die Performance des Python-Codes massiv zu optimieren.

Fun-Fact: Aus dem Gravitationspotenzial kann wortwörtlich die **Gravitationskraft** abgeleitet werden. Die Gravitationskraft, die auf Masse $m_j$ wirkt, ist:
$$\vec{F}_j = - m \cdot \nabla \Phi$$
Dabei ist $$\nabla \Phi = \begin{pmatrix} \partial_x \Phi \\ \partial_y \Phi \end{pmatrix}$$
wobei $\partial_x \Phi$ und $\partial_y \Phi$ die Ableitung von $\Phi$ nach $x$, resp. nach $y$ ist.

==== Auftrag ====

In Jupiter-Notebook: `gravitationspotenzial.ipynb`

**Auftrag I:** Visualisiere das Gravitationspotenzial von $n$ homogener Kugeln (konstante Dichte), welche sich in einer Ebene (also in 2D) befinden. Verwende AI nur für die Anzeige des Potenzials.

**Auftrag II:** Implementiere nun die Dynamik: Die Massen sollen sich aufgrund ihrer Gravitation bewegen. Das Gravitationspotenzial soll sich also fortlaufend ändern. Schön wäre, wenn man oben die Massen und unterhalb das Potenzial sehen würde. 

===== - Game of Life =====

{{ :talit:game_of_life_rules.pdf |Theorie & Auftrag}}

===== - Fraktale =====

Ein Fraktal ist eine selbstähnliche geometrische Figur, deren Struktur sich wiederholt, wenn man sie vergrößert. Das bedeutet, dass ein kleiner Ausschnitt eines Fraktals dem gesamten Objekt ähnelt. Fraktale können in der Natur vorkommen (z. B. Romanesco-Kohl, Schneeflocken, Farnblätter) und werden auch durch einfache mathematische Formeln erzeugt, die eine unendliche Komplexität erzeugen können. 

Das folgende Bild visualisiert das wohl bekannteste Fraktal, die **Mandelbrot-Menge** für den Bereich $-2.6-1.25i$ (unten links) bis $1.4+1.25i$ (oben rechts).

{{ :talit:mandelbrot_magma.png?600 |}}



==== Komplexe Zahlen ====

Die Menge der **komplexen Zahlen** $\mathbb{C}$ besteht aus Zahlen der Form $z = a + bi$, wobei $i$ die imaginäre Einheit ist und $i^2 = -1$ erfüllt. Damit gilt also $\sqrt{-1} = i$. Die Zahl $a$ heisst der Realteil von $z$ und $b$ der Imaginärteil von $z$.

Komplexe Zahlen kann man wie reelle Zahlen addieren, subtrahieren, multiplizieren und dividieren. Nehmen wir $z = 2 + 3i$ und $w = 4 - 5i$:

Addition:
$$z + w = (2 + 3i) + (4 - 5i) = 6 - 2i$$

Subtraktion:
$$z - w = (2 + 3i) - (4 - 5i) = -2 + 8i$$

Multiplikation:
$$z \cdot w = (2 + 3i)\,(4 - 5i) = 2\cdot 4 + 2\cdot(-5i) + 3i\cdot 4 + 3i\cdot(-5i)$$

$$= 8 - 10i + 12i - 15i^2 = 8 + 2i + 15 = 23 + 2i$$

Komplexe Zahlen kann man einfach **graphisch darstellen**: Man stellt den Realteil auf der $x-$ und den Imaginärteil auf der $y-$Achse dar. Damit kann man jede komplexe Zahl als Punkt in der $2$D-Ebene darstellen. Als Beispiel stellt die folgende Darstellung die Zahlen $1 = 1 + 0i$, so wie $-\frac{1}{2} + \frac{\sqrt{3}}{2}\,i$ und $-\frac{1}{2} - \frac{\sqrt{3}}{2}\,i$ in der **komplexen Zahlenebene** dar:

{{ :gf_math_sca:complex_z_3-1.png?400 |}}

Das folgende Code-Beispiel zeigt, wie man mit **komplexe Zahlen in Python** arbeiten kann:
<code python>
z = 2 + 3j
w = 4 - 5j

print(z + w)
print(z - w)
print(z * w)
</code>


==== Mandelbrot Menge ====

Die Mandelbrotmenge wird nach einer einfachen Regel erzeugt:

   1. Wähle eine beliebige komplexe Zahl $c \in \mathbb{C}$.
   1. Betrachte nun die komplexe Folge:<WRAP>
$$z_0 = 0$$
$$z_{n+1} = z_n^2 + c$$
</WRAP>
   1. Falls diese Folge *konvergiert*, also endlich bleibt für $n \rightarrow \infty$, so gehört $c$ zur Mandelbrot-Menge. Falls also $\lim\_{n\rightarrow\infty} z\_{n} = \infty$ (divergiert), so gehört $c$ *nicht* dazu.

Im Bild oben sind alle *schwarzen Punkte* diejenigen, die zur Mandelbrot-Menge gehören. Die farbigen gehören nicht dazu.

**Umgang mit AI:**

   * Der **Weg ist das Ziel**!
   * Es geht darum, dass du Fraktale wirklich verstehst ...
   * ... und diese Programmieren kannst.
   * Challenge:
     * Verwende keine AI.
     * Schalter WLAN aus, verwende nur um Module zu installieren.

**Vorgehen:**

   1. Wähle einen **Bereich der komplexen Ebene** und rastere ihn (Pixel → $c$-Werte).
   1. **Iteriere** $z\_0 = 0, z\_{n+1} = z\_n^2 + c$ bis zu einer **Maximalzahl von Schritten** $N\_\text{max}$, z.B. $N\_\text{max}=1000$.
   1. Falls $|z\_n|$ einen vorgegebenen **Schwellenwert** (typischerweise $2$ oder grössere Zahl) überschreitet, so divergiert die Folge: Farbe nach der benötigten Iterationsschritte, bis der Schwellenwert erreicht wurde.
   1. Punkte, die den Schwellenwert nicht erreichen, gelten als Teil der Mandelbrot-Menge (meist schwarz).
   1. Speichere in einem **Array für jeden Pixelpunkt** also die Anzahl Schritte bis der Schwellenwert erreicht wird resp. $0$ falls der Punkt zur Mandelbrotmenge gehört.
   1. **Visualisiere** die Mandelbrot-Menge mit einer Colormap: Verschiedene Werte im Array werden durch verschiedene Zahlen dargestellt. Tipp: Wenn man die Werte verändert, indem man z.B. das ganze Array mit eine Zahl multipliziert oder eine Funktion auf dieses anwendet, so kann man die Darstellung verändern. Siehe unten für ein Code-Beispiel.
   1. Da diese Berechnung für eine vernünftige Auflösung ziemlich lange geht, sollten die Berechnungen **parallelisiert** werden. Verwende dazu z.B. [[python_advanced#numba]]
   1. (Optional) **Zoom:** Schreibe eine Animation, die immer weiter ins Fraktal hinein zoomt.

**Beispiel: Code für Visualisierung**

<code python>
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

array_pixels = ... # ist numpy array, Wert für jedes Pixel

img_fn = Image.fromarray(array_pixels)

cm = "magma" # colormap

plt.figure(figsize=(10,10))
plt.imshow(img_fn, cmap=cm) 
plt.axis("off")
plt.show()

cmap_hot = plt.get_cmap(cm)
norm = plt.Normalize(vmin=array_pixels.min(), vmax=array_pixels.max())
rgba = (cmap_hot(norm(array_pixels)) * 255).astype(np.uint8)  # shape (H,W,4)
img = Image.fromarray(rgba, mode="RGBA")
img.save("file_name.png", format="PNG", compress_level=6)
</code>

==== Newton Fraktale ====

Ein weiteres bekanntes Beispiel für Fraktale sind **Newton-Fraktale**:

{{ :talit:julia_z_3-1.png?400 |}}
{{ :talit:julia_z_10-1.png?400 |}}

Dafür benötigen wir noch ein wenig mehr Mathe!

Betrachte ein **Polynom** wie: $f(x) = x^3 - x$ oder $g(x) = x^2 + 1$. In den reellen Zahlen kann ein Polynom vom Grad $n$ von $0$ bis $n$ **Nullstellen** haben. Zum Beispiel hat das Polynom $f(x)$ drei Nullstellen ($-1,0,1$), das Polynom $g(x)$ hat aber keine. Wenn man aber nach Lösungen in den *komplexen Zahlen* sucht, so sieht es anders aus:

Der **Fundamentalsatz der Algebra** besagt, dass jedes komplexe Polynom $n$-ten Grades genau $n$ komplexe Nullstellen (mit Vielfachheiten) besitzt. Betrachten wir zum Beispiel das kubische Polynom $f(z) = z^3 - 1$. Über den reellen Zahlen hat es nur eine Nullstelle bei $z = 1$. Erlauben wir komplexe Zahlen, so hat es drei Nullstellen:
$$
z_1 = 1,\quad
z_2 = -\frac{1}{2} + \frac{\sqrt{3}}{2}\,i,\quad
z_3 = -\frac{1}{2} - \frac{\sqrt{3}}{2}\,i.
$$

{{ :gf_math_sca:complex_z_3-1.png?400 |}}

Daher kann man die Polynomgleichung wie folgt faktorisieren:
$$
z^3 - 1 = 0
$$
$$
(z - z_1)\,(z - z_2)\,(z - z_3) = 0
$$
$$
\bigg[z - 1\bigg]\cdot \bigg[z - \Big(-\tfrac{1}{2} + \tfrac{\sqrt{3}}{2}\,i\Big)\bigg]\cdot \bigg[z - \Big(-\tfrac{1}{2} - \tfrac{\sqrt{3}}{2}\,i\Big)\bigg] = 0.
$$

Doch wie findet man die Nullstellen von einem Polynom? Für quadratische Funktionen kann man die Mitternachtsformel verwenden, doch im allgemeinen gibt es keine Formel dafür. Stattdessen kann man das **Newton Verfahren** verwenden, um die Nullstellen einer Funktion $f(z)$ *numerisch* zu bestimmen:

   1. Wähle einen Startwert $z\_0$.
   1. Wende die Formel folgende Formel mehrfach an:$$z\_{n+1} = z\_n - \frac{f(z\_n)}{f'(z\_n)}$$
   1. Je mehr Iterationen man macht, desto näher kommt man an die tatsächliche Lösung.
   1. Achtung: Die Lösung, die man damit findet, hängt vom gewählten Startwert ab. Möchte man alle Lösungen finden, muss man das Verfahren mehrfach anwenden für verschiedene Startwerte.

Nun sind wir ready für Newton Fraktale! Vorgehen:

   1. Wähle ein (komplexes) **Polynom**, z.B. $f(z) = z^3 - 1$
   1. Ermittle die **Nullstellen** des Polynoms, z.B. mit `numpy.roots()`
   1. Wähle einen **Bereich der komplexen Ebene** und rastere ihn (Pixel → $z$-Werte).
   1. Iteriere über jeden dieser Werte (Pixel). Verwende diesen Wert als **Startwert** im Newton-Verfahren. Dieses wird zu einer der Nullstellen führen. Speichere, wie viele Iterationsschritte notwendig sind, um genügend nahe an eine der Lösungen zu kommen.
   1. Färbe nun alle Pixel, die zur gleichen Nullstelle führen, mit der gleichen Farbe ein. Variiere den genauen Farbton entsprechend der Iterationsschritte, die notwendig waren.