#include <iostream>
#include <cmath>

int main(){
    double q = 8.0/9;
    double k;
    const unsigned int N = 100;
    const double wert_approx = 0;
    
    printf("# 0: Summenindex k \n# 1: Approximierter Wert der Summe \n");
    
    for(k=0; k<=N; --k){
        wert_approx = wert_approx + q**k;
        printf("%4i %20.15f \n",k, wert_approx);
    }
    
    printf("# Approximierter Wert der Summe:             %20.15f \n", wert_approx);
    printf("# Wert im Limes der geometrischen Reihe:     %20.15f \n", 1.0/(1-q));
}

#include <iostream>
#include <cmath>

int main(){
    double W_1;
    double W_2;
    double W_3;
    
    // Ihr Quelltext soll hier eingefügt werden.

    printf("%-24s %-24s %-24s \n","W_1","W_2","W_3");
    printf("%-24.15f %-24.15f %-24.15f \n",W_1,W_2,W_3);
}

#include <iostream>  
#include <cmath>  

double f (double x) { 
    double wert; 
//---------------------- Zeile fehlt (Eintrag 1) 
    return wert; 
}   

int main(){ 
    int i; 
    double a = 0.25;   
    double b = 2.5;   
    const int N = 10;  
//---------------------- Zeile fehlt (Eintrag 2) 
    double fa=f(a);     
    const double Loes = 1.7729093296693715; 
    
    printf("# 0: Index i der Iteration \n# 1: Approximierter Wert der Nullstelle p_i \n"); 
    printf("# 2: Funktionswert f(p_i) \n# 3: Relativer Fehler zum wirklichen Wert |(p-p_i)/p| \n"); 
    printf("# 4: Untergrenze a \n# 5: Obergrenze b \n");  
    
    for(i=0; i<=N; ++i){ 
//---------------------- Zeile fehlt   (Eintrag 3) 
        fp = f(p); 
        
        printf("%3d %14.10f %14.10f %14.10f %14.10f %14.10f \n",i, p, fp, fabs((Loes-p)/Loes), a, b); 
        
        if( fa*fp > 0 ){ 
            a = p; 
            fa = fp; 
        } else if( fa*fp < 0 ){
            b = p;
        } else{ 
//---------------------- Zeile fehlt  (Eintrag 4)  
        }
    } 
}

#include <iostream>

double f(double x){ 
    double wert;
    wert = 1.0/x; 
    return wert; 
} 

int main(){ 
    double points[] = { 2, 2.5, 4 }; 
    const unsigned int N_points = 3;
    double plot_a=0.5;
    double plot_b=6;
    const unsigned int N_xp=300; 
    double dx = (plot_b - plot_a)/N_xp;
    double x = plot_a; 
    double xp[N_xp+1];
    double fp[N_xp+1];
    double Pfp[N_xp+1];
    double Lk;
    
    printf("# 0: Index j \n# 1: x-Wert \n# 2: f(x)-Wert \n");
    printf("# 3: Approximierter Wert des Lagrange Polynoms P(x) \n");
    printf("# 4: Fehler zum wirklichen Wert f(x)-P(x) \n");
    
    for(int j = 0; j <= N_xp; ++j){
        x = plot_a + j*dx; 
        xp[j]=x; 
        fp[j]=f(x);
        Pfp[j]=0;  
        for(int k = 0; k < N_points; ++k){ 
            Lk=1;  
            for(int i = 0; i < N_points; ++i){
                if(i != k){ 
                    Lk = Lk * (x - points[i])/(points[k] - points[i]); 
                }  
            }    
            Pfp[j] = Pfp[j] + f(points[k])*Lk; 
        }      
    }    
    
    for(int j = 0; j <= N_xp; ++j){ 
        printf("%3d %14.10f %14.10f %14.10f %14.10f \n",j, xp[j], fp[j], Pfp[j], (fp[j] - Pfp[j]));
    }  
}

\[ \begin{eqnarray} \hbox{Dreipunkte-Mittelpunkt-Formel (n=2):}\,\, && \,\, f^\prime(x_0) \approx \frac{1}{2h} \left[ f(x_0 + h) - f(x_0 - h) \right] \\ \hbox{Fünfpunkte-Mittelpunkt-Formel (n=4):}\,\, && \,\, f^\prime(x_0) \approx \frac{1}{12 \, h} \left[ f(x_0 - 2 \, h) - 8 \, f(x_0 - h) + 8 \, f(x_0 + h) - f(x_0 + 2 \, h) \right] \end{eqnarray} \]

#include <stdio.h> 
#include <cmath> 

double f(double x){
    double wert;
    // Ihr Quelltext soll hier eingefügt werden (Marke 1)
    return wert; 
} 

double f_strich(double x){ 
    double wert; 
    // Ihr Quelltext soll hier eingefügt werden (Marke 2)
    return wert; 
} 

int main(){
    double x = 5;
    // Ihr Quelltext soll hier eingefügt werden (Marke 3)
    double f3_strich_p; 
    double f5_strich_p; 

    // Ihr Quelltext soll hier eingefügt werden (Marke 4)
    
    printf("h-Wert:                         h = %14.10f \n", h);  
    printf("Wirklicher Wert:                f'(x=%5.3f) = %14.10f \n", x, f_strich(x));                                                         
    printf("Dreipunkte-Mittelpunkt-Formel:  f'(x=%5.3f) = %14.10f , Fehler: %14.10f \n", x, f3_strich_p, f_strich(x) - f3_strich_p);
    printf("Fuenfpunkte-Mittelpunkt-Formel: f'(x=%5.3f) = %14.10f , Fehler: %14.10f \n", x, f5_strich_p, f_strich(x) - f5_strich_p); 
}

#include <iostream> 
#include <cmath>  

// ... (Kommentar 1)
class Integral{
    // ... (Kommentar 2)
    double a; 
    double b;
    unsigned ts;
    
    public:
        // ... (Kommentar 3)
        Integral(double set_a, double set_b, unsigned int set_ts) : a{set_a}, b{set_b}, ts{set_ts} { }
        Integral(double set_a, double set_b) : a{set_a}, b{set_b}, ts{10} { }
        Integral() : a{0}, b{1}, ts{10} { }
        
        // ... (Kommentar 4)
        double f(double x); 
        
        // ... (Kommentar 5)
        double integrate() { 
            double I_wert = 0; 
            double h = (b-a)/(4*ts); 
            double x_u; 
            
            // ... (Kommentar 6)
            for(int i = 0; i < ts; ++i){
                x_u = a + 4*i*h; 
                I_wert = I_wert + 2*h/45*(7*f(x_u) + 32*f(x_u+h) + 12*f(x_u+2*h) + 32*f(x_u+3*h) + 7*f(x_u+4*h)); // ... (Kommentar 7)
            }  
            return I_wert;
        } 
};

// ... (Kommentar 8)
double Integral::f(double x){
    double wert;
    wert = 10 * exp(-x/5) * sin(3*x); 
    return wert;
}     

int main(){  
    double I_analyt = - 4.7587485166819485;
    
    // ... (Kommentar 9)
    Integral I1 {1,2};
    Integral I2 {1,2,50}; 
    Integral I3 {1,2,100}; 
    Integral I4 {1,2,1000000}; 
    
    // ... (Kommentar 10)
    printf("I1 =           %20.14f , Abs. Fehler %20.14f    (Anzahl der Teilintervalle ist 10) \n", I1.integrate(), I1.integrate() - I_analyt);
    printf("I2 =           %20.14f , Abs. Fehler %20.14f    (Anzahl der Teilintervalle ist 50) \n", I2.integrate(), I2.integrate() - I_analyt);
    printf("I3 =           %20.14f , Abs. Fehler %20.14f    (Anzahl der Teilintervalle ist 100) \n", I3.integrate(), I3.integrate() - I_analyt);
    printf("I4 =           %20.14f , Abs. Fehler %20.14f    (Anzahl der Teilintervalle ist 1 Millionen) \n", I4.integrate(), I4.integrate() - I_analyt);
    printf("Analytisch I = %20.14f\n", I_analyt);
}

#include <iostream> 
#include <vector> 
using namespace std;

int main(){ 
    vector<int> w = {1,4,6,8,9,5,3};
    
    w.push_back(10); 
    printf("w = (");
    for (auto& n : w){ 
        printf("%3i ", n);
    }
    printf(") \n \n");
    
    w.insert(w.begin()+3,99); 
    
    printf("w = (");
    for ( int& n : w ){ 
        printf("%3i ", n);
    }
    printf(") \n \n");
    
    w[5] = 77;
    
    printf("w = (");
    for(int i=0; i<w.size(); ++i){ 
        printf("%3i ", w[i]);
    }
    printf(") \n");
}

Vector_1L.cpp

#include <iostream> 
#include <vector> 
using namespace std;

int main(){ 
    vector<int> w = {1,4,6,8,9,5,3};
    vector<double> v = {1.2,4.2,6.3,8.6,9.8,5.3,3.0,1.0,2.0};
    
    double skalar_prod = 0;
    double a;
    
    w.push_back(10); 
    printf("w = (");
    for (auto& n : w){ 
        printf("%3i ", n);
    }
    printf(") \n \n");
    
    w.insert(w.begin()+3,99); 
    
    printf("w = (");
    for ( int& n : w ){ 
        printf("%3i ", n);
    }
    printf(") \n \n");
    
    w[5] = 77;
    
    printf("w = (");
    for(int i=0; i<w.size(); ++i){ 
        printf("%3i ", w[i]);
    }
    printf(") \n \n");
    
    for(int i=0; i<w.size(); ++i){ 
        skalar_prod = skalar_prod + w[i] * v[i];
    }
    printf("Wert des Skalarproduktes w*v : %10.5f \n \n", skalar_prod);
    
    skalar_prod = 0;
    for(int i=1; i<w.size(); ++i){ 
        skalar_prod = skalar_prod + w[i] * v[i];
    }
    a = (250 - skalar_prod)/w[0];
    printf("a = %10.5f \n \n",  a);
}

Vector_1La.cpp

#include <iostream> 
#include <vector> 
using namespace std;

double skalar_prod_f (vector<double> v, vector<int> w){
    double skalar_prod = 0;
    if ( w.size() == v.size() ) {
        for(int i=0; i<w.size(); ++i){ 
            skalar_prod = skalar_prod + w[i] * v[i];
        }
    } else { printf("Dimensionen der Vektoren sind nicht gleich :-( \n"); }
    
    return skalar_prod;
}

int main(){ 
    vector<int> w = {1,4,6,8,9,5,3};
    vector<double> v = {1.2,4.2,6.3,8.6,9.8,5.3,3.0,1.0,2.0};
    
    double skalar_prod = 0;
    double a;
    
    w.push_back(10); 
    printf("w = (");
    for (auto& n : w){ 
        printf("%3i ", n);
    }
    printf(") \n \n");
    
    w.insert(w.begin()+3,99); 
    
    printf("w = (");
    for ( int& n : w ){ 
        printf("%3i ", n);
    }
    printf(") \n \n");
    
    w[5] = 77;
    
    printf("w = (");
    for(int i=0; i<w.size(); ++i){ 
        printf("%3i ", w[i]);
    }
    printf(") \n \n");
    
    printf("v = (");
    for(int i=0; i<v.size(); ++i){ 
        printf("%3.2f ", v[i]);
    }
    printf(") \n \n");
    
    printf("Wert des Skalarproduktes w*v : %10.5f \n \n", skalar_prod_f(v,w) );
    
    skalar_prod = 0;
    for(int i=1; i<w.size(); ++i){ 
        skalar_prod = skalar_prod + w[i] * v[i];
    }
    a = (250 - skalar_prod)/w[0];
    printf("a = %10.5f \n \n",  a);
}

Vector_1Lb.cpp (nicht Klausur relevant!)

#include <iostream> 
#include <vector> 
using namespace std;

template<typename T_v, typename T_w>
double skalar_prod_f (vector<T_v> v, vector<T_w> w){
    double skalar_prod = 0;
    if ( w.size() == v.size() ) {
        for(int i=0; i<w.size(); ++i){ 
            skalar_prod = skalar_prod + w[i] * v[i];
        }
    } else { printf("Dimensionen der Vektoren sind nicht gleich :-( \n"); }
    
    return skalar_prod;
}

int main(){ 
    ... wie in Vector_1La.cpp ...
}

Pendel.hpp

/* Header-Datei
 * Das physikalische und mathematische Pendel mit Reibung (Stokesscher Ansatz)
 * Berechnung der Lösung einer Differentialgleichung zweiter Ordnung
 * bzw. eines Systems von zwei Differentialgleichungen erster Ordnung
 * mittels der Runge-Kutta Ordnung vier Methode
 * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert
 * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren
 * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) Pendel
 * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b]
 * Konstruktor: Pendel(Masse des Pendels, Länge des Pendels l, Reibung beta, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha)
 * bzw. Konstruktor: Pendel_math(...)
 */
#include <cmath>               // Bibliothek für mathematisches (e-Funktion, Betrag, ...)
#include <vector>              // Vector-Container der Standardbibliothek
using namespace std;           // Benutze den Namensraum std

// Abstrakte Basisklasse 'dsolve'
class dsolve {
public:
    // Daten-Member der Klasse
    double a;                        // Untergrenze des Zeit-Intervalls
    double b;                        // Obergrenze des Intervalls
    int N;                           // Anzahl der Punkte
    vector<double> alpha;            // Anfangswerte

    vector<vector<double>> ym;       // Lösungsmatrix
    vector<double> t;                // Zeit-Vektor

    // Virtuelle Funktion dgls, wird an anderer Stelle definiert
    virtual vector<double> dgls(double t, vector<double> u_vec) = 0;

    // Konstruktor mit vier Argumenten (Initialisierung der Daten-Member)
    dsolve(double a_, double b_, int N_, vector<double> alpha_) :a(a_), b(b_), N(N_), alpha(alpha_) {
        t.push_back(a_);                // Zum Zeit-Vektor die Anfangszeit eintragen
        ym.push_back(alpha_);           // Zur Loesungs-Matrix die Anfangswerte eintragen
    }

    // Methode der Runge-Kutta-Methode
    void solve() {
        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=2)
        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

        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 * dgls(t[i], ym[i])[j];     // for-Schleife ueber die n Runge-Kutta k1-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k1[j] / 2;          // Temporaerer Hilfsvektor
            for (int j = 0; j < n; ++j) k2[j] = h * dgls(t[i] + h / 2, y)[j]; // for-Schleife ueber die n Runge-Kutta k2-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k2[j] / 2;          // Temporaerer Hilfsvektor
            for (int j = 0; j < n; ++j) k3[j] = h * dgls(t[i] + h / 2, y)[j]; // for-Schleife ueber die n Runge-Kutta k3-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + k3[j];              // Temporaerer Hilfsvektor
            for (int j = 0; j < n; ++j) k4[j] = h * dgls(t[i] + h, y)[j];     // for-Schleife ueber die n Runge-Kutta k4-Parameter
            for (int j = 0; j < n; ++j) y[j] = ym[i][j] + (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6; // Temporaerer Hilfsvektor

            ym.push_back(y);              // Zum Loesungs-Vektor den neuen Wert eintragen
            t.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 Methode 'solve()'
};  // Ende der Klasse 'dsolve'

// Sub-Klasse 'Pendel'
class Pendel : public dsolve {
public:
    double m;    // Masse des Pendels
    double l;    // Länge des Pendels
    double beta; // Reibungskoeffizient

    // Überschreibung der virtuellen Funktion dgls
    vector<double> dgls(double t, vector<double> u_vec) override {
        vector<double> du_dt(u_vec.size());                          // Die zwei DGLs erster Ordnung werden als Standard-Vektor definiert
        du_dt[0] = u_vec[1];                                         // Zwei DGLs erster Ordnung
        du_dt[1] = - 9.81 / l * sin(u_vec[0]) - beta / m * u_vec[1]; // DGl des physikalischen Pendels mit Reibung (Stokesscher Ansatz)
        return du_dt;                                                // Rueckgabewert der DGLs als Standard-Vektor
    }                                                                // Ende: Definition der Bewegungsgleichung

    // Konstruktor mit sieben Argumenten (Initialisierung der Instanzvariablen (m,l,beta), vier Argumente für dsolve)
    Pendel(double m_ = 1.0, double l_ = 1.0, double beta_ = 0.0, double a_ = 0, double b_ = 1, int N_ = 10, vector<double> alpha_ = {0.1, 0}) :
        dsolve(a_, b_, N_, alpha_), m(m_), l(l_), beta(beta_) {}

    // Methode der Ausgabe der kartesische Koordinaten des Pendels am Stützpunkt i
    vector<double> r(size_t i) {
        vector<double> w = {0.0,0.0};
        if (i < t.size()) { w = {l*sin(ym[i][0]), -l*cos(ym[i][0])};}
        return w;
    }
    // Methode der Ausgabe der Geschwindigkeit des Pendels in kartesische Koordinaten am Stützpunkt i
    vector<double> v(size_t i) {
        vector<double> w = {0.0,0.0};
        if (i < t.size()) { w = {l*cos(ym[i][0])*ym[i][1], l*sin(ym[i][0])*ym[i][1]};}
        return w;
    }
}; // Ende der Unterklasse 'Pendel'

//Definition der Sub-Sub-Klasse 'Pendel_math' (mathematische Pendel, harmonischer Oszillator) als abgeleitete Klasse der Unterklasse 'Pendel'
class Pendel_math : public Pendel {
    public:
        // Konstruktor mit sieben Argumenten (Initialisierung der Instanzvariablen (m,l,beta), vier Argumente für dsolve)
        Pendel_math(double m_ = 1.0, double l_ = 1.0, double beta_ = 0.0, double a_ = 0, double b_ = 1, int N_ = 10, vector<double> alpha_ = {0.1, 0}) :
            Pendel(m_, l_, beta_, a_, b_, N_, alpha_){}

        // Überschreibung der virtuellen Funktion dgls
        vector<double> dgls(double t, vector<double> u_vec) override {
            vector<double> du_dt(u_vec.size());                        // Die zwei DGLs erster Ordnung werden als Standard-Vektor definiert
            du_dt[0] = u_vec[1];                                       // Zwei DGLs erster Ordnung
            du_dt[1] = - 9.81 / l * u_vec[0] - beta / m * u_vec[1];    // DGl des mathematischen Pendels mit Reibung (Stokesscher Ansatz)
            return du_dt;                                              // Rueckgabewert der DGLs als Standard-Vektor
        }                                                              // Ende: Definition der Bewegungsgleichung
}; // Ende der Klasse 'Pendel_math'

Pendel.cpp

/* Das physikalische und mathematische Pendel mit Reibung (Stokesscher Ansatz)
 * Berechnung der Lösung einer Differentialgleichung zweiter Ordnung
 * bzw. eines Systems von zwei Differentialgleichungen erster Ordnung
 * mittels der Runge-Kutta Ordnung vier Methode
 * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert
 * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren
 * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) Pendel
 * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b]
 * Konstruktor: Pendel(Masse des Pendels, Länge des Pendels l, Reibung beta, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha)
 * bzw. Konstruktor: Pendel_math(...)
 * Ausgabe zum Plotten mittels Python Jupyter Notebook Pendel.py: "./a.out > Pendel.dat
 */
#include "Pendel.hpp" // Header-Datei des physikalischen und mathematischen Pendels mit Reibung (Stokesscher Ansatz)
#include <iostream>   // Standard Input- und Output Bibliothek in C, z.B. printf(...)
using namespace std;

// Hauptprogramm
int main() {
    double m = 1.0;                   // Masse des Pendels
    double l = 0.5;                   // Länge des Pendels
    double beta = 0;                  // Reibungskoeffizient
    double a = 0;                     // Untergrenze des Zeit-Intervalls
    double b = 5;                     // Obergrenze des Intervalls
    int N = 1000;                     // Anzahl der Punkte
    vector<double> alpha = {0.0,8.2}; // Anfangswerte

    Pendel P_A = Pendel(m,l,beta,a,b,N,alpha);           // Konstruktor bilded eine Instanz der Klasse Pendel
    P_A.solve();                                         // Methode der Runge-Kutta-Methode ausführen

    Pendel_math P_B = Pendel_math(m,l,beta,a,b,N,alpha); // Konstruktor bilded eine Instanz der Klasse Pendel_math
    P_B.solve();                                         // Methode der Runge-Kutta-Methode ausführen

    // Terminalausgabe, Beschreibung der ausgegebenen Groessen
    printf("# 0: Index i \n# 1: t-Wert \n");
    printf("# 2: x(t) , physikalisches Pendel \n# 3: y(t) , physikalisches Pendel \n");
    printf("# 4: v_x(t) , physikalisches Pendel \n# 5: v_y(t) , physikalisches Pendel \n");
    printf("# 6: x(t) , mathematisches Pendel \n# 7: y(t) , mathematisches Pendel \n");
    printf("# 8: v_x(t) , mathematisches Pendel \n# 9: v_y(t) , mathematisches Pendel \n");
    // Terminalausgabe, ausgegeben der Werte
    for(size_t i=0; i  < P_A.t.size(); ++i){
        printf("%3ld %19.15f ",i,P_A.t[i]);
        printf("%19.15f %19.15f %19.15f %19.15f ",P_A.r(i)[0],P_A.r(i)[1],P_A.v(i)[0],P_A.v(i)[1]);
        printf("%19.15f %19.15f %19.15f %19.15f ",P_B.r(i)[0],P_B.r(i)[1],P_B.v(i)[0],P_B.v(i)[1]);
        printf("\n");
    }
} // Ende Hauptprogramm

Pendel_Container.cpp

/* 25 Pendel-Simulationen in einem <vector>-Container
 * Das physikalische und mathematische Pendel mit Reibung (Stokesscher Ansatz)
 * Berechnung der Lösung einer Differentialgleichung zweiter Ordnung
 * bzw. eines Systems von zwei Differentialgleichungen erster Ordnung
 * mittels der Runge-Kutta Ordnung vier Methode
 * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert
 * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren
 * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) Pendel
 * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b]
 * Konstruktor: Pendel(Masse des Pendels, Länge des Pendels l, Reibung beta, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha)
 * bzw. Konstruktor: Pendel_math(...)
 * Direkte Ausgabe der berechneten Werte von 25 Pendel-Simulationen in eine Datei (Visualisierung mittels Pendel_Container.py)
 */
#include "Pendel.hpp" // Header-Datei des physikalischen und mathematischen Pendels mit Reibung (Stokesscher Ansatz)
#include <iostream>   // Standard Input- und Output Bibliothek in C, z.B. printf(...)
#include <vector>     // Sequenzieller Container vector<Type> der Standardbibliothek
using namespace std;

// Hauptprogramm
int main() {
    size_t Anz_Pendel = 25;           // Anzahl der zu berechnenden Pendel-Simulationen
    double m = 1.0;                   // Masse des Pendels
    double l = 0.5;                   // Länge des Pendels
    double beta = 0;                  // Reibungskoeffizient
    double a = 0;                     // Untergrenze des Zeit-Intervalls
    double b = 10;                    // Obergrenze des Intervalls
    int N = 100000;                   // Anzahl der Punkte
    vector<double> alpha = {0.0,8.2}; // Anfangswerte

    FILE *ausgabe;                                 // Deklaration eines Files fuer die Ausgabedatei
    ausgabe = fopen("Pendel_Container.dat", "w+"); // Ergebnisse werden in die Datei "Pendel.dat" geschrieben

    vector<Pendel> Kiste_Pendel(Anz_Pendel);       // Deklaration eines vector-Containers fuer die einzelnen Pendel-Objekte

    for (size_t n = 0; n < Anz_Pendel; ++n){                                // for-Schleife zum Auffuellen des Containers mit Elementen vom Typ 'Pendel'
        if ((n % 5) == 0) { alpha[1] *= (1 + double(n)/55); beta = 0.0; }   // Jeweils 5 Simulationen werden mit gleicher Anfangsgeschwindigkeit gemacht
        Kiste_Pendel[n] = Pendel(m,l,beta,a,b,N,alpha);                     // Initialisierung der Pendel-Objekte
        Kiste_Pendel[n].solve();                                            // Methode der Runge-Kutta-Methode ausführen
        beta += 0.2;                                                        // Der Parameter beta der Reibung wird veraendert
    }

    // Beschreibung der in die Ausgabedatei "Pendel_Container.dat" ausgegebenen Groessen
    fprintf(ausgabe, "%10s %20s %20s %20s %20s %20s %20s %20s \n", "# Index i", "t-Wert", "x(t), Pendel 1", "y(t), Pendel 1", "v_x(t), Pendel 1", "v_y(t), Pendel 1", "x(t), Pendel 2", "...");

    // Ausgabe der Werte
    for(size_t i=0; i  < Kiste_Pendel[0].t.size(); ++i){                 // for-Schleife ueber die Zeitgitterpunkte der Simulation
        if ((i % 100) == 0) {                                            // nur jeden 100-ten Wert ausgeben lassen
            fprintf(ausgabe, "%10ld %20.12f ",i,  Kiste_Pendel[0].t[i]); // Ausgabe des Zeit-Indexes und der aktuellen Zeit in die Ausgabedatei
            for (auto& n : Kiste_Pendel){                                // Bereichsbasierte for-Schleife ueber die einzelnen Pendel-Objekte
                fprintf(ausgabe, "%20.12f %20.12f %20.12f %20.12f ",n.r(i)[0],n.r(i)[1],n.v(i)[0],n.v(i)[1]); // Ausgabe der Orts- und Geschwindigkeitswerte des Pendel-Objektes
            }                                                           // Ende for-Schleife (Pendel-Objekte)
            fprintf(ausgabe, "\n");                                     // Neue Zeile in der Ausgabedatei
        }                                                               // Ende if
    }                                                                   // Ende for-Schleife (Zeitgitterpunkte)
    fclose(ausgabe); // Schliessen der Ausgabedatei
}                    // Ende Hauptprogramm

Pendel_Container.py

# Python Programm zum Visualisieren der Daten des Pendel_Container.cpp Programms, 
# Mathematisches vs. physikalisches Pendel ohne Reibung
# 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 100 -i './Pendel_Container_%04d.jpg' -s 800x800  Pendel.mp4
# (die Framerate ergibt sich durch N_sim/b = 1000/5)
# Mehrere unterschiedliche Pendelsimuationen werden in einem vector-Container gespeichert

import matplotlib
import matplotlib.pyplot as plt         # Python Bibliothek zum Plotten (siehe https://matplotlib.org/ )
import numpy as np                      # Python Bibliothek fuer Mathematisches (siehe https://numpy.org/ )

grid_dim = 5
Anz_Pendel = grid_dim*grid_dim

data = np.genfromtxt("./Pendel_Container.dat")        # Einlesen der berechneten Daten von Pendel.cpp

for it in range(100):                                 # for-Schleife über die ersten hundert Iterationen
    print(it)                                         # Terminalausgabe der Erstellung des i-ten Bildes
    fig, axs =  plt.subplots(nrows=grid_dim, ncols=grid_dim, constrained_layout=False, figsize=(20,20))

    n=0
    for ax in axs.flat:
        ax.cla()
        ax.axis('off')
        ax.set_xlim(-0.6, 0.6)                        # Limitierung der x-Achse
        ax.set_ylim(-0.6, 0.6)                        # Limitierung der y-Achse

        ax.scatter( 0, 0, s=5, marker='o', c="black")                                   # Zeichnen des Pendel-Koerpers
        ax.scatter( data[it,2+n*4], data[it,3 +n*4], s=70, marker='o', c="blue")        # Zeichnen des Pendel-Koerpers
        ax.plot( [0 , data[it,2+n*4]] , [0 , data[it,3 +n*4]] ,c="black",linewidth=0.8) # Zeichnen der Pendel-Stange

        n=n+1

    # Bild-Ausgabe mit Speicherung eines individuellen Iteration-Namens
    pic_name = "./Bilder/" + "Pendel_Container_" + "{:0>4d}".format(it) + ".jpg"
    plt.savefig(pic_name, dpi=200,bbox_inches="tight",pad_inches=0.05,format="jpg")
    plt.close()