Die schwingende Kette
Das Projekt Die schwingende Kette ist ein Beispiel eines schwingenden, gekoppelten Massensystems und ein zentrales Problem aus der klassischen Mechanik (siehe Walter Greiner, 'Klassische Mechanik II' [8. Auflage, 2008, Kapitel I7. Seite 76]). Das System besteht aus einem masselosen Faden, der mit $N+1$ Massenpunkten (Dingen, Teilchen, Perlen) besetzt ist. Die $N+1$ 'Perlen' sollen in diesem Projekt als Objekte einer Klasse programmiert werden. Jede Perle besitzt eine ganzzahlige Instanzvariable $n$ (die Nummer der Perle), wobei die erste ($n=0$) und letzte Perle ($n=N+1$) so befestigt werden, dass sie sich nicht im Raum bewegen können. Beschreibt man das System im zweidimensionalen Raum und spannt die erste Perle bei ($x_0=0$ ,$y_0=0$) und die letzte Perle bei ($x_{N+1}=L$ ,$y_{N+1}=0$) ein, so soll eine über den gesamten Faden konstante Fadenspannung $T$ entstehen. Die Perlen seien in äquidistanten Abständen $a$ auf dem Faden angeordnet und die Auslenkung aus der Ruhelage in y-Richtung sei relativ klein, sodass die geringfügige Auslenkung in x-Richtung zu vernachlässigbar ist. Betrachten Sie somit zunächst den Fall einer ausschließlich vertikalen Auslenkung der Teilchen (Perlen).
Die y-Auslenkung des $n$-ten Teilchens $y_n$ wird durch die Auslenkung des ($n-1$)-ten Teilchens $y_{n-1}$ und des ($n+1$)-ten Teilchens $y_{n+1}$ beeinflusst. Die rücktreibenden Kräfte $\vec{F}_{n-1}$ und $\vec{F}_{n+1}$ besitzen die folgende Form \[ \begin{eqnarray} \vec{F}_{n-1} &=& - \left( T \cdot \hbox{sin}(\alpha_{n-1}) \right) \, \vec{e}_y \underbrace{\approx}_{\hbox{lineare Näherung}} - T \left( \frac{y_n - y_{n-1}}{a} \right) \, \vec{e}_y \\ \vec{F}_{n+1} &=& - \left( T \cdot \hbox{sin}(\alpha_{n+1}) \right) \, \vec{e}_y \underbrace{\approx}_{\hbox{lineare Näherung}} - T \left( \frac{y_n - y_{n+1}}{a} \right) \, \vec{e}_y \quad , \end{eqnarray} \] wobei $\alpha_{n-1}$ der mit der Horizontalen eingeschlossene Auslenkungswinkel zum ($n-1$)-ten Teilchen und $\alpha_{n+1}$ der Winkel zu seinem anderen Nachbarn, dem($n+1$)-ten Teilchen ist. Die Bewegungsgleichung des $n$-ten Teilchens in der linearen Näherung lautet somit: \[ \begin{eqnarray} m_n \, \frac{d^2 y_n}{dt^2} \, \vec{e}_y = \vec{F}_{n-1} + \vec{F}_{n+1} &=& - T \left( \frac{y_n - y_{n-1}}{a} \right) \, \vec{e}_y - T \left( \frac{y_n - y_{n+1}}{a} \right) \, \vec{e}_y \qquad \Leftrightarrow \\ \frac{d^2 y_n(t)}{dt^2} &=& \frac{T}{m_n \, a} \left( y_{n-1} - 2 \, y_n + y_{n+1} \right) \quad , \end{eqnarray} \] wobei $m_n$ die Masse des $n$-ten Teilchens ist und der Index $n$ von $n=1$ bis $n=N$ läuft. Die Bewegungsgleichung der schwingenden Kette stellt somit ein System von $N$ Differentialgleichungen zweiter Ordnung dar. Die Bewegungsgleichung für die erste und letzte freie Perle ($n=1$ und $n=N$) vereinfacht sich, da die ganz äußeren Eckperlen eingespannt sind (Randbedingung: $y_0(t)=0$ und $y_{N+1}(t)=0 \, , \,\, \forall \, t \in \,$ℝ) \[ \begin{equation} \frac{d^2 y_1(t)}{dt^2} = \frac{T}{m_1 \, a} \left( - 2 \, y_1 + y_{2} \right) \qquad \hbox{und} \qquad \frac{d^2 y_N(t)}{dt^2} = \frac{T}{m_N \, a} \left( y_{N-1} - 2 \, y_N \right) \quad . \end{equation} \]
Analytische Lösung mit Python
Bevor wir die numerische Lösung simulieren betrachten wir uns zunächst die allgemeine Lösung des Systems von gekoppelten DGLs am Beispiel einer schwingenden Kette mit drei Perlen ($N=3$). Beim Klicken auf das nebenstehende Bild gelangen Sie zu dem Jupyter Notebook SchwingendeKette.ipynb in dem die analytischen und numerischen Lösungen in Python programmiert sind. Die analytische Lösung des Jupyter Notebooks basiert auf dem Python Modul "sympy", das ein Computer-Algebra-System für Python bereitstellt und eine Vielzahl an symbolischen Berechnungen im Bereich der Mathematik und Physik relativ einfach möglich macht.
Zunächst definieren wir uns das Systems von gekoppelten DGLs als sympy-Gleichung. Wir setzen der Einfachheit halber die Fadenspannung $T=1$ und die äquidistanten Abstände zwischen den Perlen auf $a=1$. Zusätzlich nehmen wir an, dass die Perlen alle die gleiche Masse hätten und setzen $m_1=m_2=m_3=1$. Mittels des sympy-Befehls "dsolve_system()" bestimmen wir die allgemeine analytische Lösung des Systems von gekoppelten DGLs. Die allgemeine analytische Lösung besitzt noch sechs unbestimmte Konstanten $C_1,C_2,C_3,C_4,C_5,C_6$, die durch die sechs Anfangsbedingung (Anfangsorte und Anfangsgeschwindigkeiten der drei Perlen) bestimmt sind. Wir setzen: $y_1(0)=\alpha_1 , y_2(0)=\alpha_2 , y_3(0)=\alpha_3$ , $\dot{y}_1(0)=\alpha_4 ,\dot{y}_2(0)=\alpha_5$ und $\dot{y}_3(0)=\alpha_6$ und geben diese allgemeine analytische Lösung des Systems von gekoppelten DGLs zum Vergleich mit unserem C++ Programm aus.
Um einen Einblick in die Bewegung der drei Perlen zu erhalten, müssen wir die Anfangsbedingungen festlegen. Wir lenken zum Beispiel zur Zeit $t=0$ nur die erste Perle um Eins aus ($y_1(0)=1$ , $y_2(0)=0$ , $y_3(0)=0$ , $\dot{y}_1(0)=0$ , $\dot{y}_2(0)=0$ und $\dot{y}_3(0)=0$). Die untere Abbildung veranschaulicht die Bewegung der schwingenden Kette in einer Animation. Wir erstellen diese Animation mittels des Python Moduls Matplotlib und benutzen hierbei die Funktion "FuncAnimation(...)" (siehe Animations using Matplotlib).
In dem Jupyter Notebook SchwingendeKette.ipynb berechnen wir auch zusätzlich die Eigenfrequenzen und Eigenschwingungen der schwingenden Kette auf analytischem Wege (für drei und fünf Perlen). Eine Eigenschwingung der Kette liegt dann vor, wenn alle Perlen mit der gleichen Frequenz $\omega$ schwingen. Die untere nebenstehende Abbildung veranschaulicht eine solche Eigenschwingung am Beispiel einer schwingenden Kette bestehend aus fünf Perlen.
Numerische Lösung mit C++
Umschreiben der Bewegungsgleichung in ein System von Differentialgleichungen erster Ordnung
Bevor wir uns mit dem C++ Programm zum numerischen Lösen der Bewegung der schwingenden Kette näher befassen, schreiben wir das System von $N$ Bewegungsgleichung zweiter Ordnung der schwingenden Kette in ein System von $2 N$ Differentialgleichungen erster Ordnung um. Unser System von Differentialgleichung zweiter Ordnung lautet wie folgt \begin{equation} \frac{d^2 y_n(t)}{dt^2} = \frac{T}{m_n \, a} \left( y_{n-1}(t) - 2 \, y_n(t) + y_{n+1}(t) \right) \,\, \quad \hbox{bzw.} \quad \ddot{y_n} = \frac{T}{m_n \, a} \left( y_{n-1} - 2 \, y_n + y_{n+1} \right) =: f_n(y_{n-1},y_n,y_{n+1}) \quad, \end{equation} wobei $n \in [1,N]$. Die befestigten Eckpunkte (Eckperlen) der schwingenden Kette werden durch $n=0$ und $n=N+1$ beschrieben. Wir machen die vorgegebene Variablenumbenennung ($u_n(t)=y_n(t)$, $u_{N+n}(t)=\frac{y_n}{dt}(t)=\dot{y}_n(t)$) und definieren das System von DGLs wie folgt: \begin{eqnarray} &\dot{u}_1(t) = \frac{d u_1}{dt} = \frac{d y_1}{dt} = u_{N+1}(t)& \\ &\dot{u}_2(t) = \frac{d u_2}{dt} = \frac{d y_2}{dt} = u_{N+2}(t) &\\ &... = ... &\\ &\dot{u}_N(t) = \frac{d u_N}{dt} = \frac{d y_N}{dt} = u_{2N}(t) &\\ &\dot{u}_{N+1}(t) = \frac{d u_{N+1}}{dt} = \frac{d \dot{y}_1}{dt} = \ddot{y}_1(t) =: f_1(y_{0},y_1,y_{2}) = f_1(u_{0},u_1,u_{2})&\\ &\dot{u}_{N+2}(t) = \frac{d u_{N+2}}{dt} = \frac{d \dot{y}_2}{dt} = \ddot{y}_2(t) =: f_2(y_{1},y_2,y_{3}) = f_2(u_{1},u_2,u_{3})&\\ &... = ... &\\ &\dot{u}_{2N}(t) = \frac{d u_{2N}}{dt} = \frac{d \dot{y}_N}{dt} = \ddot{y}_N(t) =: f_N(y_{N-1},y_N,y_{N+1}) = f_2(u_{N-1},u_N,u_{2N+1})& \quad , \end{eqnarray} wobei die befestigten Eckpunkte (Eckperlen) der schwingenden Kette durch $u_{0}$ und $u_{2N+1}$ beschrieben werden. Wir berechnen nun die numerische Lösung $\vec{u}(t) = \left( u_1(t), u_2(t), ..., u_{2N}(t) \right)$ in einem bestimmten Zeitintervall (z.B. $t \in [0,1]$) bei vorgegebenen Anfangsbedingungen $$ \begin{eqnarray} &u_n(0) = y_n(0) = \alpha_n \,\, , \forall n \in [1,N] & \\ &u_n(0) = \dot{y}_{n-N}(0) = \alpha_n \,\, , \forall n \in [N+1,2N] & \end{eqnarray} $$
Das C++ Programm
Im Unterpunkt Differentialgleichungen: Numerische Lösung von Anfangswertproblemen hatten wir mittels des C++ Programmes DGL_1.cpp die numerische Lösung einer Differentialgleichung erster Ordnung berechnet. Dieses Programm erweiterten wir dann im Unterpunkt Systeme von gekoppelten Differentialgleichungen und Differentialgleichungen zweiter Ordnung um auch Systemen von gekoppelten Differentialgleichungen erster Ordnung und Differentialgleichungen zweiter Ordnung numerisch Lösen zu können (siehe DGL_2_a.cpp).
Im Hinblick, auch auf die noch folgenden Projekte, ist es sicherlich hilfreich zunächst den Kern-Algorithmus der Berechnung der numerischen Lösung (Runge-Kutta Ordnung vier Verfahren) als eine C++ Klasse auslagern, die wir dann zum simulieren der schwingende Kette verwenden werden. Das folgende C++ Programm stellt z.B. eine Klasse dar, welche die Lösung der folgenden Differentialgleichung erster Ordnung $$ \begin{equation} \dot{y}(t) = \frac{d y(t)}{dt} = y - t^2 + 1 \quad, \, \hbox{mit:} \,\, a=0 \leq t \leq b=4 \, , \,\, y(a)=y(0)=\alpha=0.3 \end{equation} $$ mittels des Runge-Kutta Ordnung vier Verfahrens generiert.
// Klasse dsolve /* Berechnung der Loesung einer Differentialgleichung der Form y'=f(t,y) * mittels Runge-Kutta Ordnung vier Verfahren * Verfahren zur Loesung der DGL ist in eine Klasse ausgelagert * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] * Konstruktor: dsolve(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswert alpha=y(a)) */ #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 class dsolve{ //Definition der Klasse 'dsolve' double a = 0; // Untergrenze des Zeit-Intervalls [a,b] in dem die Loesung berechnet werden soll double b = 2; // Obergrenze des Intervalls [a,b] int N = 10; // Anzahl der Punkte in die das t-Intervall aufgeteilt wird double alpha = 0.5; // Anfangswert bei t=a: y(a)=alpha vector<double> y_RungeK_4; // Deklaration eines double Vektors zum speichern der Loesung vector<double> Zeit; // Deklaration eines double Vektors zum speichern der Zeit-Werte public: // Konstruktor mit vier Argumenten (Initialisierung der Parameter, Berechnung der Loesung der DGL) dsolve(double a_, double b_, int N_, double alpha_) : a(a_),b(b_),N(N_),alpha(alpha_) { double h = (b - a)/N; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) double k1, k2, k3, k4; // Deklaration der vier Runge-Kutta Parameter Zeit.push_back(a_); // Zum Zeit-Vektor die Endzeit eintragen y_RungeK_4.push_back(alpha_); // Zum y-Vektor den Anfangswert alpha_=y(a) eintragen for(int 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.push_back(y_RungeK_4[i] + (k1 + 2*k2 + 2*k3 + k4)/6); // Zum Loesungs-Vektor den neuen Wert eintragen Zeit.push_back(a + (i+1)*h); // Zum Zeit-Vektor die neue Zeit eintragen } // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls }; // Ende des Konstruktors double f(double t, double y); // Deklaration der Member-Funktion f(t,y) (Definition findet ausserhalb der Klasse statt) double y_analytisch(double t, double alpha); // Deklaration der Member-Funktion f(t,alpha) (Definition findet ausserhalb der Klasse statt) const vector<double>& get_y() const { return y_RungeK_4; } // Definition der konstanten Member-Funktion get_y(), Rueckgabewert vector der Loesung der DGL const vector<double>& get_zeit() const { return Zeit; } // Definition der konstanten Member-Funktion get_zeit(), Rueckgabewert vector der zeit-Punkte }; // Ende der Klasse double dsolve::f(double t, double y){ // Definition der Funktion f(t,y), Gleichung der DGL dy/dt=f(t,y) double wert; wert = y - pow(t,2) + 1; // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(t,y) } // Ende der Funktion f(t,y) double dsolve::y_analytisch(double t, double alpha){ // Analytische Loesung der DGL double wert; wert = (alpha + (pow(t,2) + 2*t + 1)*exp(-t) -1)*exp(t); // Eigentliche Definition der analytische Loesung return wert; // Rueckgabewert } // Ende der Definitiom int main(){ // Hauptfunktion dsolve Loes {0,4,30,0.3}; // Benutzt Konstruktor mit a=0, b=4, N=30 und alpha=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(int i=0; i < Loes.get_zeit().size(); ++i){ // for-Schleife zur Terminalausgabe der Loesung printf("%3d %19.15f %19.15f %19.15f %19.15e \n",i, Loes.get_zeit()[i], Loes.get_y()[i], Loes.y_analytisch(Loes.get_zeit()[i], 0.3), Loes.get_y()[i] - Loes.y_analytisch(Loes.get_zeit()[i], 0.3)); } } // Ende der Hauptfunktion
Das oben abgebildete C++ Programm zeigt zunächst eine Möglichkeit der Umformulierung des C++ Programmes DGL_1.cpp in eine Klassenstruktur, wobei als Klassenname 'dsolve' verwendet wurde. Mittels des Konstruktors werden die Daten-Member der Klasse neu initialisiert und außerdem werden die, als private Daten-Member deklarierten Standardvektoren 'vector<double> y_RungeK_4' und 'vector<double> Zeit' im Anweisungsblock des Konstruktors mit den Lösungswerten gefüllt. Der Benutzer kann dann vom Hauptprogramm mittels der public definierten Memberfunktionen 'const vector<double>& get_y()' und 'const vector<double>& get_zeit()' auf die berechneten Lösungs- und Zeitwerte zugreifen.
Nun möchten wir diese Klasse auch zum Lösen eines Systems von Differentialgleichungen benutzen. Im Unterpunkt Systeme von gekoppelten Differentialgleichungen und Differentialgleichungen zweiter Ordnung hatten wir dies für ein Systemen von zwei gekoppelten Differentialgleichungen erster Ordnung vorgestellt (siehe DGL_2.cpp). Da wir bei der schwingenden Kette ein System von $2 N$ Differentialgleichungen erster Ordnung lösen wollen, müssen wir die Klasse 'dsolve' allgemeiner gestalten und die Anzahl der zu lösenden Differentialgleichungen sollte im Programm frei wählbar sein. Das folgende C++ Programm stellt eine allgemeine Klasse 'dsolve' dar und man kann sie beliebig zum Lösen von Systemen von gekoppelten Differentialgleichungen und Differentialgleichungen zweiter Ordnung verwenden. Im unteren Programm wird sie verwendet um die Bewegung der schwingenden Kette mit drei Perlen zu simulieren (siehe SchwingendeKette_3P.cpp). Man kann sie jedoch auch für nur eine Differentialgleichung erster Ordnung benutzen (siehe dsolve.cpp) oder eine schwingende Kette mit 200 Perlen simulieren (siehe SchwingendeKette_200P.cpp).
// Die schwingende Kette mit drei Perlen /* ist ein Beispiel eines schwingenden, gekoppelten Massensystems aus dem Bereich der Mechanik * Berechnung der Lösung eines Systems von Differentialgleichungen * mittels Runge-Kutta Ordnung vier Verfahren * Verfahren zur Loesung der DGL ist in eine Klasse ausgelagert * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] * Konstruktor: dsolve(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha=y(a), DGL als Funktion) * Ausgabe zum Plotten mittels Python Jupyter Notebook SchwingendeKette.ipynb: "./a.out > SchwingendeKette_3P.dat */ #include <iostream> // 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 #include <functional> // Funktionen in der Argumentenliste von Funktionen using namespace std; // Benutze den Namensraum std class dsolve{ //Definition der Klasse 'dsolve' double a = 0; // Untergrenze des Zeit-Intervalls [a,b] in dem die Loesung berechnet werden soll double b = 1; // Obergrenze des Intervalls [a,b] int N = 10; // Anzahl der Punkte in die das t-Intervall aufgeteilt wird vector<double> alpha = {1,0,0,0,0,0}; // Sechs Anfangswerte (Orte und Geschwindigkeiten der drei Perlen) bei t=a vector<vector<double>> y_RungeK_4; // Deklaration eine double Vektor-Matrix zum speichern der Loesungen vector<double> Zeit; // Deklaration eines double Vektors zum speichern der Zeit-Werte public: // Konstruktor mit fuenf Argumenten (Initialisierung der Parameter, Uebergabe der DGL als Funktion, Berechnung der Loesung der DGL) dsolve(double a_, double b_, int N_, vector<double> alpha_, function< vector<double>(double, vector<double>)> f) : a(a_),b(b_),N(N_),alpha(alpha_) { double h = (b - a)/N; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) int n = alpha_.size(); // Anzahl der Anfangswerte bzw. DGL-Gleichungen (hier n=6) vector<double> k1(n),k2(n),k3(n),k4(n); // Deklaration der vier Runge-Kutta Parameter fuer jede Loesung vector<double> y(n); // Temporaerer Hilfsvektor Zeit.push_back(a_); // Zum Zeit-Vektor die Anfangszeit eintragen y_RungeK_4.push_back(alpha_); // Zum Loesungs-Vektor die Anfangswerte eintragen for(int i=0; i < N; ++i){ // for-Schleife ueber die einzelnen Punkte des t-Intervalls for(int j=0; j < n; ++j){k1[j] = h*f(Zeit[i],y_RungeK_4[i])[j];} // for-Schleife ueber die sechs Runge-Kutta k1-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + k1[j]/2;} // Temporaerer Hilfsvektor for(int j=0; j < n; ++j){k2[j] = h*f(Zeit[i]+h/2,y)[j];} // for-Schleife ueber die sechs Runge-Kutta k2-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + k2[j]/2;} // Temporaerer Hilfsvektor for(int j=0; j < n; ++j){k3[j] = h*f(Zeit[i]+h/2,y)[j];} // for-Schleife ueber die sechs Runge-Kutta k3-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + k3[j];} // Temporaerer Hilfsvektor for(int j=0; j < n; ++j){k4[j] = h*f(Zeit[i]+h,y)[j];} // for-Schleife ueber die sechs Runge-Kutta k4-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])/6;} // Temporaerer Hilfsvektor y_RungeK_4.push_back(y); // Zum Loesungs-Vektor den neuen Wert eintragen Zeit.push_back(a + (i+1)*h); // Zum Zeit-Vektor die neue Zeit eintragen } // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls }; // Ende des Konstruktors const vector<vector<double>>& get_y() const { return y_RungeK_4; } // Definition der konstanten Member-Funktion get_y(), Rueckgabewert Vektor-Matrix der Loesung der DGL const vector<double>& get_zeit() const { return Zeit; } // Definition der konstanten Member-Funktion get_zeit(), Rueckgabewert vector der zeit-Punkte }; // Ende der Klasse // Definition der Bewegungsgleichung der schwingenden Kette bestehend aus drei Perlen vector<double> dgls(double t, vector<double> u_vec){ vector<double> du_dt(u_vec.size()); // Die sechs DGLs erster Ordnung werden als Standard-Vektor definiert double set_T = 1; // Fadenspannung T double set_a = 1; // Aequidistante Abstaenden zwischen den Perlen double set_m_1 = 1; // Masse der ersten Perle double set_m_2 = 1; // Masse der zweiten Perle double set_m_3 = 1; // Masse der dritten Perle du_dt[0] = u_vec[3]; // Sechs DGLs erster Ordnung du_dt[1] = u_vec[4]; // ... du_dt[2] = u_vec[5]; du_dt[3] = set_T/(set_m_1*set_a)*(u_vec[1] - 2*u_vec[0]); du_dt[4] = set_T/(set_m_2*set_a)*(u_vec[0] - 2*u_vec[1] + u_vec[2]); du_dt[5] = set_T/(set_m_3*set_a)*(u_vec[1]-2*u_vec[2]); return du_dt; // Rueckgabewert der DGLs als Standard-Vektor } // Ende: Definition der Bewegungsgleichung // Analytische Loesung (siehe Python Jupyter Notebook) vector<double> y_analytisch(double t, vector<double> alpha){ // Analytische Loesung in Abhaengigkeit der sechs Anfangswerte alpha vector<double> y; // Vektor {y_1,y_2,y_3} der Loesungen y.push_back((1.0/2.0)*(alpha[0] - alpha[2])*cos(M_SQRT2*t) + (1.0/4.0)*M_SQRT2*(alpha[3] - alpha[5])*sin(M_SQRT2*t) + (1.0/4.0)*(alpha[0] - M_SQRT2*alpha[1] + alpha[2])*cos(t*sqrt(M_SQRT2 + 2)) + (1.0/4.0)*(alpha[0] + M_SQRT2*alpha[1] + alpha[2])*cos(t*sqrt(2 - M_SQRT2)) + (1.0/8.0)*sqrt(M_SQRT2 + 2)*(alpha[3]*(2 - M_SQRT2) + 2*alpha[4]*(1 - M_SQRT2) + alpha[5]*(2 - M_SQRT2))*sin(t*sqrt(M_SQRT2 + 2)) + (1.0/8.0)*sqrt(2 - M_SQRT2)*(alpha[3]*(M_SQRT2 + 2) + 2*alpha[4]*(1 + M_SQRT2) + alpha[5]*(M_SQRT2 + 2))*sin(t*sqrt(2 - M_SQRT2))); y.push_back((1.0/8.0)*(-2*pow(M_SQRT2 + 2, 2)*(-alpha[4]*(M_SQRT2 + 2) + (1 + M_SQRT2)*(alpha[3] + alpha[5]))*sin(t*sqrt(M_SQRT2 + 2)) - pow(M_SQRT2 + 2, 2)*(2*alpha[3]*(1 + M_SQRT2) - alpha[4]*(2 - M_SQRT2)*pow(M_SQRT2 + 2, 2) + 2*alpha[5]*(1 + M_SQRT2))*sin(t*sqrt(M_SQRT2 + 2)) + 2*pow(M_SQRT2 + 2, 7.0/2.0)*((-M_SQRT2*alpha[0] + 2*alpha[1] - M_SQRT2*alpha[2])*cos(t*sqrt(M_SQRT2 + 2)) + (M_SQRT2*alpha[0] + 2*alpha[1] + M_SQRT2*alpha[2])*cos(t*sqrt(2 - M_SQRT2)) + sqrt(2 - M_SQRT2)*(alpha[3]*(1 + M_SQRT2) + alpha[4]*(M_SQRT2 + 2) + alpha[5]*(1 + M_SQRT2))*sin(t*sqrt(2 - M_SQRT2))))/pow(M_SQRT2 + 2, 7.0/2.0)); y.push_back((1.0/2.0)*(-alpha[0] + alpha[2])*cos(M_SQRT2*t) - 1.0/4.0*M_SQRT2*(alpha[3] - alpha[5])*sin(M_SQRT2*t) + (1.0/4.0)*(alpha[0] - M_SQRT2*alpha[1] + alpha[2])*cos(t*sqrt(M_SQRT2 + 2)) + (1.0/4.0)*(alpha[0] + M_SQRT2*alpha[1] + alpha[2])*cos(t*sqrt(2 - M_SQRT2)) + (1.0/8.0)*sqrt(M_SQRT2 + 2)*(alpha[3]*(2 - M_SQRT2) + 2*alpha[4]*(1 - M_SQRT2) + alpha[5]*(2 - M_SQRT2))*sin(t*sqrt(M_SQRT2 + 2)) + (1.0/8.0)*sqrt(2 - M_SQRT2)*(alpha[3]*(M_SQRT2 + 2) + 2*alpha[4]*(1 + M_SQRT2) + alpha[5]*(M_SQRT2 + 2))*sin(t*sqrt(2 - M_SQRT2))); return y; // Rueckgabe des Vektors {y_1,y_2,y_3} der Loesungen } // Ende der analytische Loesung // Hauptfunktion int main(){ vector<double> u_init {1,0,0,0,0,0}; // Vektor der sechs Anfangswerte (Orte und Geschwindigkeiten der drei Perlen) bei t=a dsolve Loes1 {0,20,100000,u_init,dgls}; // Benutzt Konstruktor mit a=0, b=20, N=100000 printf("# 0: Index i \n# 1: t-Wert \n# 2:y_1 \n# 3:y_2 \n# 4:y_3 \n"); // Beschreibung der ausgegebenen Groessen printf("# 5: Analytische Loesung y_1 \n# 6: Analytische Loesung y_2 \n# 7: Analytische Loesung y_3 \n"); // Beschreibung der ausgegebenen Groessen for(int i=0; i < Loes1.get_zeit().size(); ++i){ // for-Schleife zur Terminalausgabe der Loesung printf("%3d %19.15f %19.15f %19.15f %19.15f",i, Loes1.get_zeit()[i], Loes1.get_y()[i][0], Loes1.get_y()[i][1], Loes1.get_y()[i][2]); printf("%19.15f %19.15f %19.15f \n",y_analytisch(Loes1.get_zeit()[i],u_init)[0],y_analytisch(Loes1.get_zeit()[i],u_init)[1],y_analytisch(Loes1.get_zeit()[i],u_init)[2]); } // Ende der for-Schleife } // Ende der Hauptfunktion
Im Vergleich zu unserer ersten Klasse 'dsolve' (siehe dsolve_a.cpp) ist das zugrundeliegende, zu lösende Differentialgleichungssystem in einer vektorwertigen Funktion 'vector<double> dgls(double t, vector<double> u_vec){' definiert, die an den Konstruktor der Klasse als Argument übergeben wird. Hierfür ist die Einbindung '#include <functional>' am Anfang des Programms nötig. Mittels des Konstruktors werden die Daten-Member der Klasse neu initialisiert, wobei die Anfangswerte auch als ein Standardvektor an den Konstruktor übergeben werden (hier bei drei Perlen sechs Anfangswerte). In Abhängigkeit von der Anzahl der zu lösenden Differentialgleichungen wird die Vektor-Matrix 'vector<vector<double>> y_RungeK_4;' mit den Lösungen aufgefüllt.
Obwohl der Kern-Algorithmus zur Berechnung der numerischen Lösung von $N$-Differentialgleichungen schon allgemein in der Klasse 'dsolve' implementiert ist, kann das obere Programm nur den Fall für $N=3$ berechnen und mit der analytischen Lösung vergleichen. Bei dem folgenden Programm ist das behoben und mittels dieses C++ Programms kann man die schwingende Kette mit beliebig vielen Perlen simulieren.
// Die schwingende Kette mit beliebiger Anzahl von Perlen /* ist ein Beispiel eines schwingenden, gekoppelten Massensystems aus dem Bereich der Mechanik * Berechnung der Lösung eines Systems von Differentialgleichungen * mittels Runge-Kutta Ordnung vier Verfahren * Verfahren zur Loesung der DGL ist in eine Klasse ausgelagert * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b] * Konstruktor: dsolve(Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha=y(a), DGL als Funktion) * Ausgabe zum Plotten mittels Python Jupyter Notebook SchwingendeKette.ipynb: "./a.out > SchwingendeKette.dat */ #include <iostream> // 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 #include <functional> // Funktionen in der Argumentenliste von Funktionen using namespace std; // Benutze den Namensraum std class dsolve{ //Definition der Klasse 'dsolve' double a = 0; // Untergrenze des Zeit-Intervalls [a,b] in dem die Loesung berechnet werden soll double b = 1; // Obergrenze des Intervalls [a,b] int N = 10; // Anzahl der Punkte in die das t-Intervall aufgeteilt wird vector<double> alpha = {1,0,0,0,0,0}; // Sechs Anfangswerte (Orte und Geschwindigkeiten; hier drei Perlen) bei t=a vector<vector<double>> y_RungeK_4; // Deklaration eine double Vektor-Matrix zum speichern der Loesungen vector<double> Zeit; // Deklaration eines double Vektors zum speichern der Zeit-Werte public: // Konstruktor mit fuenf Argumenten (Initialisierung der Parameter, Uebergabe der DGL als Funktion, Berechnung der Loesung der DGL) dsolve(double a_, double b_, int N_, vector<double> alpha_, function< vector<double>(double, vector<double>)> f) : a(a_),b(b_),N(N_),alpha(alpha_) { double h = (b - a)/N; // Abstand dt zwischen den aequidistanten Punkten des t-Intervalls (h=dt) int n = alpha_.size(); // Anzahl der Anfangswerte bzw. DGL-Gleichungen vector<double> k1(n),k2(n),k3(n),k4(n); // Deklaration der vier Runge-Kutta Parameter fuer jede Loesung vector<double> y(n); // Temporaerer Hilfsvektor Zeit.push_back(a_); // Zum Zeit-Vektor die Anfangszeit eintragen y_RungeK_4.push_back(alpha_); // Zum Loesungs-Vektor die Anfangswerte eintragen for(int i=0; i < N; ++i){ // for-Schleife ueber die einzelnen Punkte des t-Intervalls for(int j=0; j < n; ++j){k1[j] = h*f(Zeit[i],y_RungeK_4[i])[j];} // for-Schleife ueber die sechs Runge-Kutta k1-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + k1[j]/2;} // Temporaerer Hilfsvektor for(int j=0; j < n; ++j){k2[j] = h*f(Zeit[i]+h/2,y)[j];} // for-Schleife ueber die sechs Runge-Kutta k2-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + k2[j]/2;} // Temporaerer Hilfsvektor for(int j=0; j < n; ++j){k3[j] = h*f(Zeit[i]+h/2,y)[j];} // for-Schleife ueber die sechs Runge-Kutta k3-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + k3[j];} // Temporaerer Hilfsvektor for(int j=0; j < n; ++j){k4[j] = h*f(Zeit[i]+h,y)[j];} // for-Schleife ueber die sechs Runge-Kutta k4-Parameter for(int j=0; j < n; ++j){y[j] = y_RungeK_4[i][j] + (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])/6;} // Temporaerer Hilfsvektor y_RungeK_4.push_back(y); // Zum Loesungs-Vektor den neuen Wert eintragen Zeit.push_back(a + (i+1)*h); // Zum Zeit-Vektor die neue Zeit eintragen } // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls }; // Ende des Konstruktors const vector<vector<double>>& get_y() const { return y_RungeK_4; } // Definition der konstanten Member-Funktion get_y(), Rueckgabewert Vektor-Matrix der Loesung der DGL const vector<double>& get_zeit() const { return Zeit; } // Definition der konstanten Member-Funktion get_zeit(), Rueckgabewert vector der zeit-Punkte }; // Ende der Klasse // Definition der Bewegungsgleichung der schwingenden Kette bestehend aus anz_p Perlen vector<double> dgls(double t, vector<double> u_vec){ int anz_p = static_cast<int>(u_vec.size()/2); // Anzahl der Perlen vector<double> du_dt; // Die 2*anz_p DGLs erster Ordnung werden als Standard-Vektor definiert double set_T = 1; // Fadenspannung T double set_a = 1; // Aequidistante Abstaenden zwischen den Perlen double set_m = 1; // Masse der ersten Perle for(int i=0; i < anz_p; ++i){ // Erster Block von anz_p DGLs erster Ordnung du_dt.push_back(u_vec[anz_p+i]); // ... } // ... du_dt.push_back(set_T/(set_m*set_a)*(u_vec[1] - 2*u_vec[0])); // Zweiter Block von anz_p DGLs erster Ordnung for(int i=anz_p+1; i < 2*anz_p-1; ++i){ // ... mit Bewegungsgleichungen du_dt.push_back(set_T/(set_m*set_a)*(u_vec[i-anz_p-1] - 2*u_vec[i-anz_p] + u_vec[i-anz_p+1])); // ... } // ... du_dt.push_back(set_T/(set_m*set_a)*(u_vec[anz_p-2]-2*u_vec[anz_p-1])); // ... return du_dt; // Rueckgabewert der DGLs als Standard-Vektor } // Ende: Definition der Bewegungsgleichung-Funktion // Hauptfunktion (hier fuer fuenf Perlen) int main(){ vector<double> u_init {1,0,0,0,0,0,0,0,0,0}; // Vektor der 10 Anfangswerte (Orte und Geschwindigkeiten der fuenf Perlen) bei t=a int anz_p = static_cast<int>(u_init.size()/2); // Anzahl der Perlen dsolve Loes1 {0,20,1000,u_init,dgls}; // Benutzt Konstruktor mit a=0, b=20, N=1000 printf("# 0: Index i \n# 1: t-Wert \n# 2:y_1 \n# 3:y_2 \n# 4:y_3 \n# 5:y_4 \n# 6:y_5 \n"); // Beschreibung der ausgegebenen Groessen for(int i=0; i < Loes1.get_zeit().size(); ++i){ // for-Schleife zur Terminalausgabe; Zeitwerte printf("%3d %19.15f ",i, Loes1.get_zeit()[i]); // Ausgabe Index und Zeitwert for(int j=0; j < anz_p; ++j){ // for-Schleife zur Terminalausgabe; n-Perlen printf("%19.15f ",Loes1.get_y()[i][j]); // Ausgabe y-Wert der Perle j } // Ende der for-Schleife; n-Perlen printf("\n"); // Neue Zeile } // Ende der for-Schleife } // Ende der Hauptfunktion
In der Hauptfunktion des Programms definiert man den Standardvektor der Anfangsbedingungen und somit legt man auch fest aus wievielen Perlen die Kette besteht (hier aus fünf Perlen). Die Bewegungsgleichung der schwingenden Kette ist allgemein für $N$-Perlen in der Vektor-Funktion 'dgls' definiert. In der Hauptfunktion führt man die Simulation durch den Aufruf des Konstruktors aus und die Simulationsergebnisse werden dann für alle Perlen mittels 'printf()' ausgegeben.
Die untere Animation zeigt eine Simulation mit $N=200$ Perlen. Die Simulationsdaten wurden mit dem C++ Programm SchwingendeKette_200P.cpp berechnet und die Visualisierung der Resultate wurde mittels des Jupyter Notebooks SchwingendeKette.ipynb angefertigt. Man kann auch die Daten direkt mit dem folgenden Python-Skript visualisieren: PythonPlot_Kette.py.
In diesem Projekt könnte man z.B. die folgenden weiteren Aufgaben betrachten:
Betrachten Sie den Fall ungleicher Massen an einem Beispiel.
Betrachten Sie den allgemeinen Fall $\vec{F}_{n-1} = - \left( T \cdot \hbox{sin}(\alpha_{n-1}) \right) \, \vec{e}_y$ und $\vec{F}_{n+1} = - \left( T \cdot \hbox{sin}(\alpha_{n+1}) \right) \, \vec{e}_y$ und beziehen zusätzlich auch die Bewegungen der Perlen in x-Richtung in Ihre Berechnungen ein. Erstellen Sie dann ein Beispiel-Simulation mit großen Perlenausschlägen und vergleichen die Ergebnisse mit den Ergebnissen der linearen Näherung. Bauen Sie auch diesen Anwendungsfall sinnvoll in die Klassenstruktur Ihres Programmes ein.