Die Musterlösung benutzt das 'Runge-Kutta Ordnung vier'-Verfahren zur Bestimmung der numerischen Lösung der Differentialgleichung. Das folgende C++ Programm zeigt die Umformulierung des C++ Programmes DGL_1.cpp mittels einer C++-Funktion, wobei als Funktionsname 'dsolve(...)' verwendet wurde.
// Musterlösung der Aufgabe 1 des Übungsblattes Nr.9 /* Berechnung der Lösung einer Differentialgleichung der Form y'=f(t,y) * mittels Runge-Kutta Ordnung vier Verfahren * Verfahren zur Loesung der DGL ist in eine Funktion dsolve(a,b,N,Anfangswert) ausgelagert * Rückgabetyp vector<vector<double>>: Matrix der y_t-Werte * Zeitentwicklung fuer aequidistante t-Werte in [a,b] */ #include <stdio.h> // Standard Input- und Output Bibliothek in C, z.B. printf(...) #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) #include <vector> // Vector-Container der Standardbibliothek using namespace std; // Benutze den Namensraum std double f(double t, double y){ // Definition der Funktion f(t,y), Gleichung der DGL dy/dt=f(t,y) double wert; wert = y - t*t + 1; // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(t,y) } // Ende der Funktion f(t,y) double y_analytisch(double t, double alpha){ // Analytische Loesung der DGL double wert; wert = t*t + 2*t + 1 + (alpha - 1) * exp(t); // Eigentliche Definition der analytische Loesung return wert; // Rueckgabewert } // Ende der Definitiom // Definition der Funktion 'dsolve' vector<vector<double>> dsolve( double a, double b, unsigned N, double alpha){ double h = (b - a)/N; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) vector<double> y_RungeK_4( N + 1 ); // Deklaration eines double Vektors zum speichern der y-Werte vector<double> Zeit( N + 1 ); // Deklaration eines double Vektors zum speichern der Zeit-Werte vector<vector<double>> y_zeit; // Deklaration der Vektor-Matrix des Rückgabe-Objektes der Funktion (y-Zeit-Werte) double k1, k2, k3, k4; // Deklaration der vier Runge-Kutta Parameter Zeit[0] = a; // Initialisierung des Zeit-Anfangswertes y_RungeK_4[0] = alpha; // Initialisierung des Anfangswertes bei t=a: y(a)=alpha for (unsigned i = 0; i < N; ++i){ // for-Schleife ueber die einzelnen Punkte des t-Intervalls k1 = h*f(Zeit[i],y_RungeK_4[i]); // Runge-Kutta Parameter 1 k2 = h*f(Zeit[i]+h/2,y_RungeK_4[i] + k1/2); // Runge-Kutta Parameter 2 k3 = h*f(Zeit[i]+h/2,y_RungeK_4[i] + k2/2); // Runge-Kutta Parameter 3 k4 = h*f(Zeit[i]+h,y_RungeK_4[i] + k3); // Runge-Kutta Parameter 4 y_RungeK_4[i+1] = y_RungeK_4[i] + (k1 + 2*k2 + 2*k3 + k4)/6; // Zum y-Vektor den neuen Wert eintragen Zeit[i+1] = a + (i+1)*h; // Zum Zeit-Vektor den neuen Wert eintragen } // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls y_zeit.push_back(y_RungeK_4); // Einfügen der berechneten y-Werte in die Rueckgabe-Matrix y_zeit.push_back(Zeit); // Einfügen der Zeit-Werte in die Rueckgabe-Matrix return y_zeit; // Rueckgabe-Matrix der y-Zeit-Werte der Loesung der DGL } // Ende der Funktion dsolve int main(){ // Definition der Matrix der Lösungswerte und Aufruf der Funktion dsolve(a,b,N,alpha) // wobei [a,b]=[0,4]: Zeit-Intervalls, N=500: Anzahl der Punkte im t-Intervall, alpha=0.3: Anfangswert bei t=a: y(a)=alpha vector<vector<double>> y_zeit = dsolve(0,4,500,0.3); printf("# 0: Index i \n# 1: t-Wert \n# 2:Runge-Kutta Ordnung vier \n"); // Beschreibung der ausgegebenen Groessen printf("# 3: Analytische Loesung \n# 4: Fehler \n"); // Beschreibung der ausgegebenen Groessen // for-Schleife zur Terminalausgabe der Loesung for(size_t i = 0; i < y_zeit[0].size(); ++i){ printf("%3ld %19.15f %19.15f %19.15f %19.15e \n",i, y_zeit[1][i], y_zeit[0][i], y_analytisch(y_zeit[1][i], 0.3), y_zeit[0][i] - y_analytisch(y_zeit[1][i], 0.3)); } }
Der Rückgabewert der Funktion ist vom Typ her eine geschachtelte vector-Struktur (Standardvektor-Matrix: vector<vector<double>>), welche die berechneten Funktions- und Zeitwerte ($(y_i,t_i) \,\, \forall i \in [0,N]$) zurückgibt. Die Argumentenliste der Funktion 'vector<vector<double>> dsolve( double a, double b, unsigned N, double alpha){ ... }' enthält die Unter- und Obergrenze des Zeit-Intervalls [a,b], die Anzahl N der Punkte im Zeit-Intervall und als letztes Argument den Parameter $\alpha$ des Anfangswertes von y ($y(a)=y(0)=\alpha$).
In der Hauptfunktion des Programms wird die Funktion dann einfach bei der Definition der Lösungsmatrix verwendet, indem sie diese mit den berechneten Werten initialisiert: vector<vector<double>> y_zeit = dsolve(0,4,500,0.3); Die berechneten Werte werden danach im Terminal ausgegeben und mit den Werten der analytischen Lösung verglichen.
Führt man das Programm aus, so erhält man die folgende Terminalausgabe, wenn man $N=50$ Zeit-Gitterpunkte $t_i, i=0,1,2, ..., N$ benutzt.
Für $N=500$ erhält man die folgenden Werte
, und für $N=5000$ ergibt sich die folgende Terminalausgabe:
Man erkennt deutlich, dass die approximierten Werte des Lösungsvektors $y_i$ sich bei ansteigender Anzahl der Zeit-Gitterpunkte $N$ dem wirklichen Wert der analytischen Lösung immer stärker annähern. Außerdem ist ein generelles Ansteigen des Fehlers ${\cal F} = y_i - y_{analytisch}(t_i)$ im Laufe der Zeit zu beobachten. Für $t=1$ und $N=50$ ($N=5000$) erhält man beispielsweise einen Fehler ${\cal F} \approx -8.4 \cdot 10^{-7}$ (${\cal F} = \approx -8.4 \cdot 10^{-15} $), wohingegen man am Ende der Simulation bei $t=4$ den Fehler ${\cal F} \approx 4.9 \cdot 10^{-6}$ (${\cal F} = \approx 1.6 \cdot 10^{-13}$) erhält.
Man hätte die Struktur der Funktion 'dsolve(...)' auch in anderer Weise konzipieren können. Die in dem folgenden C++ Programm definierte Funktion besitzt einen anderen Aufbau und der Rückgabetyp und die Argumentenliste der Funktion unterscheiden sich von der im vorigen Programm.
// Musterlösung der Aufgabe 1 des Übungsblattes Nr.9 /* Berechnung der Lösung einer Differentialgleichung der Form y'=f(t,y) * mittels Runge-Kutta Ordnung vier Verfahren * Verfahren zur Loesung der DGL ist in eine Funktion dsolve(Zeitvektor,Anfangswert) ausgelagert * Rückgabetyp vector<double>: Vektor der berechneten y-Werte * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] */ #include <stdio.h> // Standard Input- und Output Bibliothek in C, z.B. printf(...) #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) #include <vector> // Vector-Container der Standardbibliothek using namespace std; // Benutze den Namensraum std double f(double t, double y){ // Definition der Funktion f(t,y), Gleichung der DGL dy/dt=f(t,y) double wert; wert = y - t*t + 1; // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(t,y) } // Ende der Funktion f(t,y) double y_analytisch(double t, double alpha){ // Analytische Loesung der DGL double wert; wert = t*t + 2*t + 1 + (alpha - 1) * exp(t); // Eigentliche Definition der analytische Loesung return wert; // Rueckgabewert } // Ende der Definitiom vector<double> dsolve( vector<double> Zeit, double alpha){ // Definition der Funktion 'dsolve' vector<double> y_RungeK_4( Zeit.size()); // Deklaration eines double Vektors zum speichern der Loesung double k1, k2, k3, k4; // Deklaration der vier Runge-Kutta Parameter double h = Zeit[1] - Zeit[0]; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) y_RungeK_4[0] = alpha; // Initialisierung des Anfangswertes bei t=a: y(a)=alpha for (size_t i = 0; i < Zeit.size()-1; ++i){ // for-Schleife ueber die einzelnen Punkte des t-Intervalls k1 = h*f(Zeit[i],y_RungeK_4[i]); // Runge-Kutta Parameter 1 k2 = h*f(Zeit[i]+h/2,y_RungeK_4[i] + k1/2); // Runge-Kutta Parameter 2 k3 = h*f(Zeit[i]+h/2,y_RungeK_4[i] + k2/2); // Runge-Kutta Parameter 3 k4 = h*f(Zeit[i]+h,y_RungeK_4[i] + k3); // Runge-Kutta Parameter 4 y_RungeK_4[i+1] = y_RungeK_4[i] + (k1 + 2*k2 + 2*k3 + k4)/6; // Zum y-Vektor den neuen Wert eintragen } // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls return y_RungeK_4; // Rueckgabewert vector der Loesung der DGL } // Ende der Funktion dsolve int main(){ // Hauptfunktion double a = 0; // Untergrenze des Zeit-Intervalls [a,b] in dem die Loesung berechnet werden soll double b = 4; // Obergrenze des Intervalls [a,b] int N = 500; // Anzahl der Punkte in die das t-Intervall aufgeteilt wird double h = (b - a)/N; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) double alpha = 0.3; // Anfangswert bei t=a: y(a)=alpha vector<double> Zeit(N+1); // Deklaration des double-vector-Containers mit N+1 Einträgen (Gitterpunkte des Zeitintervalls) for(size_t i = 0; i < Zeit.size(); ++i){ // for-Schleife zum Festlegen der einzelnen Elemente des Zeitvektors Zeit[i] = a + i*h; // Wert-Zuweisung an das i-te Elementes des Zeitvektors } // Ende for-Schleife vector<double> y_RungeK_4 = dsolve(Zeit,alpha); // Definition des Loesungsvektors und Aufruf der Funktion dsolve(...) printf("# 0: Index i \n# 1: t-Wert \n# 2:Runge-Kutta Ordnung vier \n"); // Beschreibung der ausgegebenen Groessen printf("# 3: Analytische Loesung \n# 4: Fehler \n"); // Beschreibung der ausgegebenen Groessen for(size_t i = 0; i < Zeit.size(); ++i){ // for-Schleife zur Terminalausgabe der Loesung printf("%3ld %19.15f %19.15f %19.15f %19.15e \n",i, Zeit[i], y_RungeK_4[i], y_analytisch(Zeit[i], 0.3), y_RungeK_4[i] - y_analytisch(Zeit[i], 0.3)); } } // Ende der Hauptfunktion
Der Rückgabewert der Funktion ist vom Typ double-Standardvektor ( vector<double> dsolve(vector<double> Zeit, double alpha){ ... } ) und die in diesem Container gespeicherten Gleitkommazahlen stellen die Elemente des Lösungsvektor $y_i$ der 'Runge-Kutta Ordnung vier'-Methode dar. Die Argumentenliste der Funktion enthält als ersten Eintrag den double-Standardvektor 'Zeit' ( vector<double> Zeit ) der Zeitgitterpunkte und das zweite Argument ist der Parameter $\alpha$ des Anfangswertes von y ($y(a)=y(0)=\alpha$). Das Hauptprogramm ist nun ein wenig umfangreicher, da zunächst der Vektor der Zeitgitterpunkte erzeugt werden muss. Die berechneten Werte unterscheiden sich natürlich in den beiden Programmversionen nicht voneinander und man erhält auch mit dem zweiten C++ Programm A9_1_func_a.cpp die oben abgebildeten Terminalausgaben.
Die folgende, dritte Variante des Programms (A9_1_func_b.cpp), definiert eine Funktion 'dsolve(...)' ( vector<double> dsolve(double a, double b, unsigned N, double alpha) { ... } ), die denselben Rückgabetyp wie die vorige Funktion besitzt, jedoch die Argumentenliste der ersten Variante verwendet. Die zugehörigen Zeitwerte werden in dieser Variante direkt bei der Ausgabe im Hauptprogramm berechnet.
/* Optimierte Musterlösung der Aufgabe 1 des Übungsblatts Nr. 9 * Berechnung der Lösung einer Differentialgleichung der Form y'=f(t,y) * mittels Runge-Kutta-Verfahren 4. Ordnung. * Die Funktion dsolve(a, b, N, alpha) gibt einen Vektor der y-Werte zurück. * Zeitentwicklung für äquidistante t-Werte in [a, b]. */ #include <iostream> // Moderne C++-Ein-/Ausgabe #include <iomanip> // Für std::fixed und std::setprecision #include <vector> // Vector-Container der Standardbibliothek #include <cmath> // Bibliothek für mathematische Funktionen (exp, pow) using namespace std; // Benutze den Namensraum std // Definition der Funktion f(t, y) für die DGL y' = f(t, y) double f(double t, double y) { return y - t*t + 1; } // Analytische Lösung der DGL double y_analytisch(double t, double alpha) { return t*t + 2*t + 1 + (alpha - 1) * exp(t); } // Definition der Funktion dsolve vector<double> dsolve(double a, double b, unsigned N, double alpha) { double h = (b - a) / N; // Schrittweite h = dt vector<double> y_RungeK_4(N + 1); // Vektor für y-Werte y_RungeK_4[0] = alpha; // Anfangswert y(a) = alpha for (unsigned i = 0; i < N; ++i) { double t = a + i * h; // Zeitpunkt t_i double k1 = h * f(t, y_RungeK_4[i]); double k2 = h * f(t + h/2, y_RungeK_4[i] + k1/2); double k3 = h * f(t + h/2, y_RungeK_4[i] + k2/2); double k4 = h * f(t + h, y_RungeK_4[i] + k3); y_RungeK_4[i + 1] = y_RungeK_4[i] + (k1 + 2*k2 + 2*k3 + k4) / 6; } return y_RungeK_4; // Rückgabe der y-Werte } int main() { double a = 0.0; // Untergrenze des Zeitintervalls double b = 4.0; // Obergrenze des Zeitintervalls unsigned N = 500; // Anzahl der Zeitschritte double alpha = 0.3; // Anfangswert y(0) = alpha double h = (b - a) / N; // Schrittweite h = dt // Numerische Lösung berechnen vector<double> y_RungeK_4 = dsolve(a, b, N, alpha); cout << fixed << setprecision(15); // 15 Nachkommastellen für Präzision cout << "# 0: Index i\n# 1: t-Wert\n# 2: Runge-Kutta Ordnung vier\n"; cout << "# 3: Analytische Lösung\n# 4: Fehler\n"; for (unsigned i = 0; i <= N; ++i) { double t = a + i * h; cout << i << " " << t << " " << y_RungeK_4[i] << " " << y_analytisch(t,alpha) << " " << scientific << y_RungeK_4[i] - y_analytisch(t,alpha) << fixed << "\n"; } }
Das folgende Programm stellt eine Abänderung des C++ Programms Vector_Dinge.cpp dar. In der inline Funktion 'Gehe_Zeitschritt(...)' wurden einerseits periodische Randbedingungen anstelle der Reflexion an den Rändern der Kiste eingebaut, andererseits wurde eine zusätzlich Teilchen/Umgebungseigenschaft in die Funktion implementiert. Es handelt sich dabei um eine fiktive Box im unteren linken Bereich der Kiste, die die Eigenschaft besitzt, dass falls ein Teilchen in diesen Bereich gerät, es ruhend, entsprechend seiner Nummer auf der x-Achse positioniert wird.
// Musterlösung der Aufgabe 2 des Übungsblattes Nr.9 /* Ausgabe zum Plotten (Python) mittels: "./a.out > A9_2.dat" * python3 A9_2.py */ #include <iostream> // Ein- und Ausgabebibliothek #include <vector> // Sequenzieller Container vector<Type> der Standardbibliothek #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) using namespace std; // Benutze den Namensraum std //Definition der Klasse 'Ding' class Ding{ // Oeffentliche Klasse public: // Instanzvariablen (Daten-Member) der Klasse unsigned int n; // Nummer des Dinges double x, y, z; // Ort des Dinges double v_x, v_y, v_z; // Geschwindigkeit des Dinges // Konstruktor mit sieben Argumenten und Standardinitialisierungen Ding(unsigned int set_n = 0, double set_x = 0.0, double set_y = 0.0, double set_z = 0.0, double set_v_x = 0.0, double set_v_y = 0.0, double set_v_z = 0.0) : n{set_n}, x{set_x}, y{set_y}, z{set_z}, v_x{set_v_x}, v_y{set_v_y}, v_z{set_v_z} {} // Member-Funktionen der Klasse void Gehe_Zeitschritt(double dt, double max_x, double max_y, double max_z); // Deklaration einer Member-Funktion (Definition findet ausserhalb der Klasse statt) }; /* Definition der Funktion Gehe_Zeitschritt(...) als inline-Methode der Klasse Ding * Die Funktion beschreibt die Veränderung der Ortskoordinaten (x,y,z) bei einem Zeitschritt dt * Es sind zusaetzlich spezielle Randbedingungen an die Bewegung formuliert, so dass sich die * Dinge nur in einem Bereich von [0,max_x], [0,max_y] und [0,max_z] bewegen koennen * periodische Randbedingungen + zusätzliche Eigenschaft in der linken unteren Ecke */ inline void Ding::Gehe_Zeitschritt(double dt, double max_x, double max_y, double max_z){ x += v_x * dt; y += v_y * dt; z += v_z * dt; // Periodische Randbedingungen der Kiste if (x < 0){x = max_x;} if (x > max_x){x = 0;} if (y < 0){y = max_y;} if (y > max_y){y = 0;} if (z < 0){z = max_z;} if (z > max_z){z = 0;} // Falls das Ding in die untere linke Ecke kommt, wird es // auf der Horizontalen bei y=0 entsprechend seiner Nummer angeordnet if (x < 5 && y < 5){ v_x = 0; v_y = 0; v_z = 0; x = max_x*n/70; y = 0; z = 0; } } int main(){ // Hauptfunktion double kiste_x = 40; // Laenge der Kiste double kiste_y = 40; // Breite der Kiste unsigned int Anz_Teilchen = 70; // Definition der Anzahl der zu erzeugenden Dinge double t; // Deklaration des Zeitparameters double dt = 0.2; // Definition der Laenge des Zeitschrittes double Anz_tSchritte = 300; // Anzahl der Zeitschritte vector<Ding> Kiste_Teilchen; // Deklaration eines vector-Containers for (unsigned int n = 0; n < Anz_Teilchen; ++n){ // for-Schleife zum Auffuellen des Containers mit Elementen vom Typ 'Ding' // Initialisierung diaginal in der x-y-Ebene, Geschwindigkeit in x- und y-Richtung Kiste_Teilchen.push_back( Ding {n, n*kiste_x/Anz_Teilchen + 0.1, n*kiste_x/Anz_Teilchen + 0.1, 0, sin(n+1.0), cos(n+1.0), 0} ); } printf("# 0: Index i \n# 1: Zeit t \n# 2: Nummer des Teilchens 1 \n"); // Beschreibung der ausgegebenen Groessen printf("# 3: x-Position des Teilchens 1 \n# 4: y-Position des Teilchens 1 \n"); // Beschreibung der ausgegebenen Groessen printf("# 5: Nummer des Teilchens 2 \n# 6: .... bis Anzahl der Teichen \n"); // Beschreibung der ausgegebenen Groessen for (int i = 0; i < Anz_tSchritte; ++i){ // for-Schleife fuer die zeitliche Entwicklung der Dinge in der Kiste t = i * dt; // Zeit geht in dt-Schritten voran printf("%3i %10.6f ", i, t); // Ausgabe Index i und Zeit t for (auto& j : Kiste_Teilchen){ // Bereichsbasierte for-Schleife zum Ausgeben der x-y-Orte der Dinge printf("%3i %10.6f %10.6f ", j.n, j.x, j.y); // Ausgabe Teilchenorte j.Gehe_Zeitschritt(dt,kiste_x, kiste_y, 0.0); // Aufruf der inline Funktion Gehe_Zeitschritt(...) } printf("\n"); } }
Mittels des Konstruktors wurden 70 Objekte der Klasse Ding mit unterschiedlichen Anfangsorten und Geschwindigkeiten erstellt. Einige der Teilchen mit niedriger Nummer $n$ sind bereits bei der Initialisierung im unteren Bereich der Kiste, sodass sie gleich beim ersten Zeitschritt auf der x-Achse ruhend positioniert werden. Das folgende Python-Skript visualisiert dann die vom C++ Programm berechneten Teilchenbewegungen.
# Python Plot zur Musterlösung der Aufgabe 2 des Übungsblattes Nr.9 # Python Programm zum Plotten der Daten des Vector-Containers mit 70 Dingen (A9_2.cpp) # Es werden hier mehrere Bilder der zeitlichen Entwicklung des Systems in einem Ordner 'Bilder' gespeichert # !!!! Sie muessen vor der Ausfuehrung des Programms den Ordner Bilder erstellen !!!! # Die einzelnen Bilder kann mann dann mittels des folgenden Terminalbefehls zu einem Video binden: # ffmpeg -framerate 5 -i './Vector_Dinge_%03d.png' -c:v libx264 A9_2.mp4 import matplotlib.pyplot as plt # Python Bibliothek zum Plotten (siehe https://matplotlib.org/ ) import matplotlib.patches as patches import numpy as np # Python Bibliothek fuer Mathematisches (siehe https://numpy.org/ ) data = np.genfromtxt("./A9_2.dat") # Einlesen der berechneten Daten von A9_2.cpp fig, ax = plt.subplots() ax.set_title(r'Container mit Teilchen') # Titel der Abbildung ax.set_ylabel('y') # Beschriftung y-Achse ax.set_xlabel('x') # Beschriftung x-Achse r = 200 # Radius eines Dings plot_min=0 # Festlegung der x-Untergrenze (Abmessung Kiste) plot_max=40 # Festlegung der x-Obergrenze (Abmessung Kiste) anz_teilchen = int(len(data[0,:])/3) # Definition der Anzahl der Dinge cmap = plt.cm.nipy_spectral # Definition der Farbschattierung der Dinge line_colors = cmap (np.linspace(0.2,1,anz_teilchen)) # Definition der Farbschattierung der Dinge alp = 0.9 # Transparenz eines Dinges Box = patches.Rectangle((0, 0), 40, 40, linewidth=3, edgecolor='grey', facecolor='none') Box_klein = patches.Rectangle((0, 0), 5, 5, linewidth=1, edgecolor='red', facecolor='grey', alpha=0.6) for it in range(len(data[:,0])): # for-Schleife fuer die zeitliche Entwicklung der Dinge in der Kiste print(it) # Terminalausgabe der Erstellung des i-ten Bildes ax.cla() ax.add_patch(Box) ax.add_patch(Box_klein) for i in range(anz_teilchen): # for-Schleife ueber die Teilchen in der Kiste ax.scatter(data[it,3*i+3],data[it,3*i+4], marker='o', color=line_colors[i], s=r) # Kennzeichnung der Position des Dinges durch einen blauen Kreis ax.text(data[it,3*i+3],data[it,3*i+4], str(int(data[it,3*i+2])), fontsize=10, verticalalignment='center', horizontalalignment='center', color="white", alpha=alp) # Ding Nr. ax.set_xlim(plot_min-1,plot_max+1) # Plot-Limit x-Achse ax.set_ylim(plot_min-1,plot_max+1) # Plot-Limit y-Achse # Bild-Ausgabe mit Speicherung eines individuellen Iteration-Namens pic_name = "./Bilder/" + "Vector_Dinge_" + "{:0>3d}".format(it) + ".png" plt.savefig(pic_name, dpi=200,bbox_inches="tight",pad_inches=0.05,format="png")
Die unten abgebildete Animation zeigt die Bewegung der 70 Teilchen und wurde mittels des oberen Python-Skriptes erstellt. Im Laufe der Simulation (die aus 300 Einzelbildern besteht) gelangen weitere Teichen in den unteren linken Bereich, und auch sie werden dann bei $y=0$ positioniert.