Systeme von gekoppelten Differentialgleichungen und Differentialgleichungen zweiter Ordnung

Im vorigen Unterpunkt hatten wir die unterschiedlichen Verfahren zum Lösen von Differentialgleichungen erster Ordnung kennengelernt. Die Bewegungsgleichungen vieler physikalischer Systeme sind jedoch von zweiter Ordnung in der Zeit und dieses Unterkapitel der Vorlesung 9 befasst sich damit, wie man solche Differentialgleichungen höherer Ordnung numerisch mittels eines C++ Programmes löst. Um ein Differentialgleichung zweiter Ordnung mittels des Computers lösen zu können, schreibt man zunächst die DGL in ein System von zwei gekoppelten Differentialgleichungen ersten Ordnung um und diese löst man dann mit den Verfahren, die in der vorigen Vorlesung behandelt wurden. In diesem Unterpunkt werden wir uns zunächst mit Systemen von gekoppelten Differentialgleichungen erster Ordnung befassen und dann das numerische Lösen von Differentialgleichungen zweiter Ordnung vorstellen.

Systeme von gekoppelten Differentialgleichungen

Wir betrachten zunächst das numerische Lösen eines Systems von $m$-gekoppelten Differentialgleichungen (DGLs) erster Ordnung der Form $$ \begin{eqnarray} \dot{y_1}(t) &=& \frac{d y_1}{dt} = f_1(t,y_1,y_2,...,y_m) \\ \dot{y_2}(t) &=& \frac{d y_2}{dt} = f_2(t,y_1,y_2,...,y_m) \\ \dot{y_3}(t) &=& ... \,\, = \\ ... &=& ... \\ \dot{y}_m(t) &=& \frac{d y_m}{dt} = f_m(t,y_1,y_2,...,y_m) \quad , \end{eqnarray} $$ wobei die zeitliche Entwicklung der Vektorfunktion $\vec{y}(t) = \left( y_1(t), y_2(t), ..., y_m(t) \right)$ in den Grenzen $a \leq t \leq b$ gesucht wird. Die $m$-Funktionen $f_i(t,y_1,y_2,...,y_m)\, , \,\, i \in [1,2,...,m]$ bestimmen das System der DGLs und somit das Verhalten der gesuchten Funktion $\vec{y}(t)$. Es wird hierbei vorausgesetzt, dass die Funktionen $f_i(t,y_1,y_2,...,y_m)$ auf einer Teilmenge ${\cal D}$ ( $\mathbb{R}^{m+1} \supseteq {\cal D}$ ) kontinuierlich definiert sind und das so definierte Anfangswertproblem "well-posed" ist und eine eindeutige Lösung $\vec{y}(t)$ existiert. Bei gegebener Anfangskonfiguration $$ \begin{equation} y_1(a) = \alpha_1 \, , \,\, y_2(a) = \alpha_2 \, , \,\, ... \,\, ,\, y_m(a) = \alpha_m \end{equation} $$ ist es dann numerisch möglich das System von gekoppelten DGLs zu lösen.

Beispiel: Numerische Lösung eines Systems von zwei gekoppelten Differentialgleichungen erster Ordnung

Wir betrachten speziell das folgende System bestehend aus zwei gekoppelten DGLs ($m=2$): $$ \begin{eqnarray} \dot{y_1}(t) &=& \frac{d y_1}{dt} = \, 3 y_1 + 2 y_2 - \left( 2 t^2 + 1 \right) \cdot e^{2t} \, =: \, f_1(t,y_1,y_2) \\ \dot{y_2}(t) &=& \frac{d y_2}{dt} = \, 4 y_1 + y_2 + \left( t^2 + 2 t - 4 \right) \cdot e^{2t} \, =: \, f_2(t,y_1,y_2) \quad , \end{eqnarray} $$ und sind an der numerischen Lösung $\vec{y}(t) = \left( y_1(t), y_2(t) \right)$ im Zeitintervall $t \in [0,1]$ interessiert. Die Anfangsbedingungen lauten $$ \begin{equation} y_1(0) = \alpha_1 = 1 \, , \,\, y_2(0) = \alpha_2 = 1 \quad . \end{equation} $$

Das Lösen dieses Systems von DGLs ist auf gleichem Wege möglich, wie man einzelne Differentialgleichungen numerisch approximiert.

Das folgende C++ Programm löst das obere System von gekoppelten DGLs mittels der Euler- und Runge-Kutta Ordnung vier Methode. Zunächst werden die beiden System bestimmenden DGL-Funktionen $f_1(t,y_1,y_2)$ und $f_2(t,y_1,y_2)$ als C++ Funktionen definiert ('double f_1(double t, double y_1, double y_2){ ...}' und 'double f_2(double t, double y_1, double y_2){ ...}'). Danach folgen die beiden analytischen Lösungen 'y_1_analytisch' und 'y_2_analytisch', welche mittels des Jupyter Notebooks DGL_2.ipynb berechnet wurden (näheres siehe weiter unten).

DGL_2.cpp
/* Berechnung der Lösung eines Systems von Differentialgleichung 
 * der Form y1'=f1(t,y1,y2) , y2'=f2(t,y1,y2)
 * mittels der Euler Methode und Runge-Kutta Ordnung vier Methode 
 * Zeitentwicklung von y1(t), y2(t) der fuer unterschiedliche t-Werte in [a,b] 
 *  Ausgabe zum Plotten mittels Python Jupyter Notebook DGL_2.ipynb: "./a.out > DGL_2.dat
 */
#include <stdio.h>             // Standard Input- und Output Bibliothek in C, z.B. printf(...)
#include <cmath>               // Bibliothek für mathematisches (e-Funktion, Betrag, ...)

double f_1(double t, double y_1, double y_2){                                // Deklaration und Definition der Funktion f_1(t,y1,y2) 
    double wert;
    wert = 3*y_1 + 2*y_2 - (2*pow(t,2) + 1)*exp(2*t);                        // Eigentliche Definition der Funktion
    return wert;                                                             // Rueckgabewert der Funktion f_1
}                                                                            // Ende der Funktion f_1

double f_2(double t, double y_1, double y_2){                                // Deklaration und Definition der Funktion f_2(t,y1,y2) 
    double wert;
    wert = 4*y_1 + y_2 + (pow(t,2) +2*t - 4)*exp(2*t);                       // Eigentliche Definition der Funktion
    return wert;                                                             // Rueckgabewert der Funktion f_2
}                                                                            // Ende der Funktion f_2

double y_1_analytisch(double t){                                             // Analytische Loesung y_1(t) 
    double wert;                                                             // bei gegebenem Anfangswert y_1(a)=1, y_2(a)=1
    wert = exp(5*t)/3 - exp(-t)/3 + exp(2*t);                                // Eigentliche Definition der analytische Loesung
    return wert;                                                             // Rueckgabewert 
}                                                                            // Ende der Definition

double y_2_analytisch(double t){                                             // Analytische Loesung y_2(t) 
    double wert;                                                             // bei gegebenem Anfangswert y_1(a)=1, y_2(a)=1
    wert = exp(5*t)/3 + 2*exp(-t)/3 + pow(t,2)*exp(2*t);                     // Eigentliche Definition der analytische Loesung
    return wert;                                                             // Rueckgabewert 
}                                                                            // Ende der Definitiom

int main(){                                                                  // Hauptfunktion
    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 = 100;                                                             // 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_1 = 1;                                                      // 1.Anfangswert bei t=a: y_1(a)=alpha_1
    double alpha_2 = 1;                                                      // 2.Anfangswert bei t=a: y_2(a)=alpha_2
    double t;                                                                // Aktueller Zeitwert
    double y_Euler_1 = alpha_1;                                              // Deklaration und Initialisierung der numerischen Loesung der Euler Methode fuer y_1
    double y_RungeK_4_1 = alpha_1;                                           // Deklaration und Initialisierung der numerischen Loesung der Runge-Kutta Ordnung vier Methode
    double k1_1,k2_1,k3_1,k4_1;                                              // Deklaration der vier Runge-Kutta Parameter fuer y_1
    double y_Euler_2 = alpha_2;                                              // Deklaration und Initialisierung der numerischen Loesung der Euler Methode fuer y_2
    double y_RungeK_4_2 = alpha_2;                                           // Deklaration und Initialisierung der numerischen Loesung der Runge-Kutta Ordnung vier Methode
    double k1_2,k2_2,k3_2,k4_2;                                              // Deklaration der vier Runge-Kutta Parameter fuer y_2
    double tmp;                                                              // Variable zum Zwischenspeichern von Ergebnissen
    
    printf("# 0: Index i \n# 1: t-Wert \n# 2: Euler Methode y1 \n# 3: Euler Methode y2 \n");   // Beschreibung der ausgegebenen Groessen
    printf("# 4: Runge-Kutta Ordnung vier y1 \n# 5: Runge-Kutta Ordnung vier y2 \n");          // Beschreibung der ausgegebenen Groessen
    printf("# 6: Analytische Loesung y1 \n# 7: Analytische Loesung y2 \n");                    // Beschreibung der ausgegebenen Groessen
    
    for(int i=0; i <= N; ++i){                                               // for-Schleife ueber die einzelnen Punkte des t-Intervalls
        t = a + i*h;                                                         // Zeit-Parameter wird um h erhoeht
        printf("%3d %19.15f %19.15f %19.15f %19.15f ",i, t, y_Euler_1, y_Euler_2, y_RungeK_4_1);   // Ausgaben der Loesungen
        printf(" %19.15f %19.15f %19.15f \n", y_RungeK_4_2, y_1_analytisch(t), y_2_analytisch(t)); // Ausgaben der Loesungen
        
        tmp = y_Euler_1 + h*f_1(t,y_Euler_1,y_Euler_2);                      // Euler Methode
        y_Euler_2 = y_Euler_2 + h*f_2(t,y_Euler_1,y_Euler_2);                // y_2 Euler Methode
        y_Euler_1 = tmp;                                                     // y_1 Euler Methode
        
        k1_1 = h*f_1(t,y_RungeK_4_1,y_RungeK_4_2);                           // Runge-Kutta Parameter k1 fuer y_1
        k1_2 = h*f_2(t,y_RungeK_4_1,y_RungeK_4_2);                           // Runge-Kutta Parameter k1 fuer y_2
        k2_1 = h*f_1(t+h/2,y_RungeK_4_1+k1_1/2,y_RungeK_4_2+k1_2/2);         // Runge-Kutta Parameter k2 fuer y_1
        k2_2 = h*f_2(t+h/2,y_RungeK_4_1+k1_1/2,y_RungeK_4_2+k1_2/2);         // Runge-Kutta Parameter k2 fuer y_2
        k3_1 = h*f_1(t+h/2,y_RungeK_4_1+k2_1/2,y_RungeK_4_2+k2_2/2);         // Runge-Kutta Parameter k3 fuer y_1
        k3_2 = h*f_2(t+h/2,y_RungeK_4_1+k2_1/2,y_RungeK_4_2+k2_2/2);         // Runge-Kutta Parameter k3 fuer y_2
        k4_1 = h*f_1(t+h,y_RungeK_4_1+k3_1,y_RungeK_4_2+k3_2);               // Runge-Kutta Parameter k4 fuer y_1
        k4_2 = h*f_2(t+h,y_RungeK_4_1+k3_1,y_RungeK_4_2+k3_2);               // Runge-Kutta Parameter k4 fuer y_2
        y_RungeK_4_1 = y_RungeK_4_1 + (k1_1 + 2*k2_1 + 2*k3_1 + k4_1)/6;     // Runge-Kutta Ordnung vier Methode fuer y_1
        y_RungeK_4_2 = y_RungeK_4_2 + (k1_2 + 2*k2_2 + 2*k3_2 + k4_2)/6;     // Runge-Kutta Ordnung vier Methode fuer y_2
    }                                                                        // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls
}                                                                            // Ende der Hauptfunktion

Der weitere Verlauf des Quelltextes ähnelt dem C++ Programm DGL_2.cpp aus dem Unterpunkt Differentialgleichungen: Numerische Lösung von Anfangswertproblemen der vorigen Vorlesung, wobei jedoch für jede der beiden Differentialgleichungen getrennt voneinander die jeweiligen Datentypen und Variablen definiert wurden. Im Terminal gibt das Programm dann die beiden gesuchten Funktionen $y_1(t)$ und $y_2(t)$ sowohl für das Euler Verfahren, als auch für die Runge-Kutta Ordnung vier Methode aus und vergleicht die berechneten Ergebnisse mit dem wirklichen Wert der analytischen Lösung.

Wieder erkennt man, dass die Abweichungen zum wirklichen Wert, am Anfang der Simulation noch recht klein sind und im Laufe der Simulation ansteigen. Auch ist die weit bessere Performance der Runge-Kutta Ordnung vier Methode gegenüber dem Euler Verfahren deutlich zu erkennen.

Differentialgleichungen zweiter Ordnung

Wir betrachten nun das numerische Lösen einer Differentialgleichung höherer Ordnung (zunächst allgemein der Ordnung $m$) mit folgender Form $$ \begin{equation} y^{(m)}(t) = \frac{d^m y(t)}{dt^m} = f(t,y,\dot{y},\ddot{y},...,y^{(m-1)})\,\, , \quad a \leq t \leq b \quad , \end{equation} $$ bei gegebener Anfangskonfiguration $$ \begin{equation} y(a) = \alpha_1 \, , \,\, \dot{y}(a) = \alpha_2 \, , \,\, ... \,\, ,\, y^{(m-1)}(a) = \alpha_m \quad . \end{equation} $$ Man kann nun diese DGL von Ordnung $m$ in ein System von $m$ Differentialgleichungen umschreiben und dieses dann numerisch lösen. Wir machen zunächst die folgende Variablenumbenennung ($u_1(t)=y(t)$, $u_2(t)=\dot{y}(t)$, $u_3(t)=\ddot{y}(t)$, ... $u_m(t)=y^{(m-1)}(t)$) und definieren das System von DGLs wie folgt: $$ \begin{eqnarray} \dot{u_1}(t) &=& \frac{d u_1}{dt} = \frac{d y}{dt} = u_2(t) \\ \dot{u_2}(t) &=& \frac{d u_2}{dt} = \frac{d \dot{y}}{dt} = \ddot{y}(t) = u_3(t) \\ \dot{u_3}(t) &=& \frac{d u_3}{dt} = \frac{d \ddot{y}}{dt} = ... \,\,\\ ... &=& ... &\\ \dot{u}_{m-1}(t) &=& \frac{d y^{(m-2)}}{dt} = y^{(m-1)}(t) = u_m(t) \\ \dot{u}_{m}(t) &=& \frac{d y_{(m-1)}}{dt} = y^{(m)}(t) = f(t,y,\dot{y},\ddot{y},...,y^{(m-1)}) = f(t,u_1,u_2,u_3,...,u_m) \quad , \end{eqnarray} $$ wobei die zeitliche Entwicklung der Vektorfunktion $\vec{u}(t) = \left( u_1(t), u_2(t), ..., u_m(t) \right)$ in den Grenzen $a \leq t \leq b$ gesucht wird. Es wird hierbei wieder vorausgesetzt, dass die Funktion $f(t,y,\dot{y},\ddot{y},...,y^{(m-1)})$ auf einer Teilmenge ${\cal D}$ ( $\mathbb{R}^{m+1} \supseteq {\cal D}$ ) kontinuierlich definiert ist und das so definierte Anfangswertproblem "well-posed" ist und eine eindeutige Lösung $\vec{u}(t)$ existiert. Bei gegebener Anfangskonfiguration $$ \begin{equation} u_1(a) = y(a) = \alpha_1 \, , \,\, u_2(a) = \dot{y}(a) = \alpha_2 \, , \,\, ... \,\, ,\, u_m(a) = y^{(m-1)}(a)= \alpha_m \end{equation} $$ ist es dann numerisch möglich, das System von gekoppelten DGLs zu lösen und somit das zeitliche Verhalten der Funktion $y(t) = u_1(t)$ zu simulieren.

Beispiel: Numerische Lösung einer Differentialgleichung zweiter Ordnung ($m=2$)

Wir betrachten speziell die folgende Differentialgleichung zweiter Ordnung ($m=2$): $$ \begin{equation} \ddot{y} - 2\dot{y} + 2y = e^{2t} \hbox{sin}(t) \,\, \quad \hbox{bzw.} \quad \ddot{y} = 2\dot{y} - 2y + e^{2t} \hbox{sin}(t) =: f(t,y,\dot{y}) \end{equation} $$ Wir machen die vorgegebene Variablenumbenennung ($u_1(t)=y(t)$, $u_2(t)=\dot{y}(t)$) und definieren das System von DGLs wie folgt: $$ \begin{eqnarray} \dot{u_1}(t) &=& \frac{d u_1}{dt} = \frac{d y}{dt} = u_2(t) \\ \dot{u_2}(t) &=& \frac{d u_2}{dt} = \frac{d \dot{y}}{dt} = \ddot{y}(t) = 2 u_2(t) - 2 u_1(t) + e^{2t} \hbox{sin}(t) =: f(t,u_1,u_2) \end{eqnarray} $$ Wir berechnen nun die numerische Lösung $\vec{u}(t) = \left( u_1(t), u_2(t) \right)$ im Zeitintervall $t \in [0,1]$ bei vorgegebener Anfangsbedingung $$ \begin{equation} u_1(0) = y(0) = \alpha_1 = -0.4 \,\, , \quad u_2(0) = \dot{y}(0) = \alpha_2 = -0.6 \quad . \end{equation} $$

Das folgende C++ Programm löst die obere Differentialgleichung zweiter Ordnung mittels der Euler- und Runge-Kutta Ordnung vier Methode.

DGL_2_a.cpp
/* Berechnung der Lösung einer Differentialgleichung zweiter Ordnung ( y''= 2*y' - 2*y + exp(2*t)*sin(t) )
 * welches in ein System von zwei Differentialgleichung umgeschrieben wurde
 * Form des DGL Systems: u_1'=y'=f_1(t,u_1,u_2)=u_2 , u_2'=y''=f_2(t,u_1,u_2) 
 * Loesung mittels der Euler Methode und Runge-Kutta Ordnung vier Methode 
 * Zeitentwicklung von u_1(t), u_2(t) fuer unterschiedliche t-Werte in [a,b] 
 *  Ausgabe zum Plotten mittels Python Jupyter Notebook DGL_2.ipynb: "./a.out >> DGL_2_a.dat
 */
#include <stdio.h>             // Standard Input- und Output Bibliothek in C, z.B. printf(...)
#include <cmath>               // Bibliothek für mathematisches (e-Funktion, Betrag, ...)

double f_1(double t, double u_1, double u_2){                                // Deklaration und Definition der Funktion f_1(t,u_1,u_2) 
    double wert;
    wert = u_2;                                                              // Eigentliche Definition der Funktion
    return wert;                                                             // Rueckgabewert der Funktion f_1
}                                                                            // Ende der Funktion f_1

double f_2(double t, double u_1, double u_2){                                // Deklaration und Definition der Funktion f_2(t,u_1,u_2) 
    double wert;
    wert = 2*u_2 - 2*u_1 + exp(2*t)*sin(t);                                  // Eigentliche Definition der Funktion
    return wert;                                                             // Rueckgabewert der Funktion f_2
}                                                                            // Ende der Funktion f_2

double u_1_analytisch(double t, double alpha_1, double alpha_2){             // Analytische Loesung y(t)=u_1(t) 
    double wert;                                                             // bei gegebenem Anfangswert y(a)=1, y'(a)=1
    wert = exp(t)*sin(t)*(alpha_2 + 1.0/5.0 - alpha_1) + exp(t)*cos(t)*(alpha_1 + 2.0/5.0) + ((sin(t) - 2*cos(t))*exp(2*t))/5.0;
    return wert;                                                             // Rueckgabewert 
}                                                                            // Ende der Definition

double u_2_analytisch(double t, double alpha_1, double alpha_2){             // Analytische Loesung y'(t)=u_2(t)
    double wert;                                                             // bei gegebenem Anfangswert y(a)=1, y'(a)=1
    wert = ((-3*cos(t) + 4*sin(t))*exp(2*t))/5.0 - 2*((alpha_1 - alpha_2/2.0 + 1.0/10.0)*sin(t) - cos(t)*(alpha_2 + 3.0/5.0)/2.0)*exp(t);
    return wert;                                                             // Rueckgabewert 
}                                                                            // Ende der Definition

int main(){                                                                  // Hauptfunktion
    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 = 100;                                                             // 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_1 = -0.4;                                                   // 1.Anfangswert bei t=a: u_1(a)=alpha_1
    double alpha_2 = -0.6;                                                   // 2.Anfangswert bei t=a: u_2(a)=alpha_2
    double t;                                                                // Aktueller Zeitwert
    double u_Euler_1 = alpha_1;                                              // Deklaration und Initialisierung der numerischen Loesung der Euler Methode fuer u_1
    double u_RungeK_4_1  =alpha_1;                                           // Deklaration und Initialisierung der numerischen Loesung der Runge-Kutta Ordnung vier Methode
    double k1_1,k2_1,k3_1,k4_1;                                              // Deklaration der vier Runge-Kutta Parameter fuer u_1
    double u_Euler_2 = alpha_2;                                              // Deklaration und Initialisierung der numerischen Loesung der Euler Methode fuer u_2
    double u_RungeK_4_2 = alpha_2;                                           // Deklaration und Initialisierung der numerischen Loesung der Runge-Kutta Ordnung vier Methode
    double k1_2,k2_2,k3_2,k4_2;                                              // Deklaration der vier Runge-Kutta Parameter fuer u_2
    double tmp;                                                              // Variable zum Zwischenspeichern von Ergebnissen
    
    printf("# 0: Index i \n# 1: t-Wert \n# 2: Euler Methode y=u1 \n# 3: Euler Methode y'=u2 \n# 4: Runge-Kutta Ordnung vier y=u1 \n"); 
    printf("# 5: Runge-Kutta Ordnung vier y'=u2 \n# 6: Analytische Loesung y=u1 \n# 7: Analytische Loesung y'=u2 \n"); // Beschreibung der ausgegebenen Groessen
    
    for(int i=0; i <= N; ++i){                                               // for-Schleife ueber die einzelnen Punkte des t-Intervalls
        t = a + i*h;                                                         // Zeit-Parameter wird um h erhoeht
        printf("%3d %19.15f %19.15f %19.15f %19.15f %19.15f ",i, t, u_Euler_1, u_Euler_2,u_RungeK_4_1, u_RungeK_4_2); // Ausgaben der Loesungen
        printf(" %19.15f %19.15f \n",u_1_analytisch(t,alpha_1,alpha_2),u_2_analytisch(t,alpha_1,alpha_2));            // Ausgaben der Loesungen

        tmp = u_Euler_1 + h*f_1(t,u_Euler_1,u_Euler_2);                      // Euler Methode
        u_Euler_2 = u_Euler_2 + h*f_2(t,u_Euler_1,u_Euler_2);                // u_2   
        u_Euler_1 = tmp;                                                     // u_1 
        
        k1_1 = h*f_1(t,u_RungeK_4_1,u_RungeK_4_2);                           // Runge-Kutta Parameter k1 fuer u_1
        k1_2 = h*f_2(t,u_RungeK_4_1,u_RungeK_4_2);                           // Runge-Kutta Parameter k1 fuer u_2
        k2_1 = h*f_1(t+h/2,u_RungeK_4_1+k1_1/2,u_RungeK_4_2+k1_2/2);         // Runge-Kutta Parameter k2 fuer u_1
        k2_2 = h*f_2(t+h/2,u_RungeK_4_1+k1_1/2,u_RungeK_4_2+k1_2/2);         // Runge-Kutta Parameter k2 fuer u_2
        k3_1 = h*f_1(t+h/2,u_RungeK_4_1+k2_1/2,u_RungeK_4_2+k2_2/2);         // Runge-Kutta Parameter k3 fuer u_1
        k3_2 = h*f_2(t+h/2,u_RungeK_4_1+k2_1/2,u_RungeK_4_2+k2_2/2);         // Runge-Kutta Parameter k3 fuer u_2
        k4_1 = h*f_1(t+h,u_RungeK_4_1+k3_1,u_RungeK_4_2+k3_2);               // Runge-Kutta Parameter k4 fuer u_1
        k4_2 = h*f_2(t+h,u_RungeK_4_1+k3_1,u_RungeK_4_2+k3_2);               // Runge-Kutta Parameter k4 fuer u_2 
        u_RungeK_4_1 = u_RungeK_4_1 + (k1_1 + 2*k2_1 + 2*k3_1 + k4_1)/6;     // Runge-Kutta Ordnung vier Methode fuer u_1
        u_RungeK_4_2 = u_RungeK_4_2 + (k1_2 + 2*k2_2 + 2*k3_2 + k4_2)/6;     // Runge-Kutta Ordnung vier Methode fuer u_2
    }                                                                        // Ende for-Schleife ueber die einzelnen Punkte des t-Intervalls
}                                                                            // Ende der Hauptfunktion

Im Terminal gibt das Programm dann die berechneten Daten für $u_1(t) = y(t)$ und $u_2(t) = \dot{y}(t)$, sowohl für das Euler Verfahren, als auch für die Runge-Kutta Ordnung vier Methode aus und vergleicht die berechneten Ergebnisse mit dem wirklichen Wert der analytischen Lösung.

Wieder erkennt man, dass die Abweichungen zum wirklichen Wert, am Anfang der Simulation noch recht klein sind und im Laufe der Simulation ansteigen. Wieder ist auch die weit bessere Performance der Runge-Kutta Ordnung vier Methode gegenüber dem Euler Verfahren deutlich zu erkennen.

Visualisierung und numerische Lösung mittels Python

Beim Klicken auf das untere rechte Bild gelangen Sie zu dem Jupyter Notebook DGL_2.ipynb in dem eine Visualisierung der Ergebnisse des oberen C++ Programms programmiert ist. Die untere Abbildung auf der linken Seite zeigt z.B. die Visualisierung der Daten von DGL_2.cpp. Zusätzlich wird in dem Notebook auch die numerische Lösung direkt in Python generiert (Methode "integrate.odeint()" im Python-Modul "scipy") und mit den simulierten Daten der unterschiedlichen Verfahren verglichen.