Das unten abgebildete C++ Programm, berechnet mittels einer modifizierten Newton-Raphson Methode die Nullstelle der Funktion $f(x) = e^x - 20$). Für den approximativen Wert der Ableitung wurde einerseits die 'Dreipunkte-Mittelpunkt-Formel' und die 'Fünfpunkte-Mittelpunkt-Formel' verwendet und zusätzlich die Ergebnisse mit einer Berechnung verglichen, die den analytischen Ausdruck für $f^\prime(x)$ benutzt.
// Musterlösung der Aufgabe 1 des Übungsblattes Nr.7 #include <iostream> // Ein- und Ausgabebibliothek #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) double f(double x) { // Deklaration und Definition der Funktion f(x) double wert; // Lokale double-Variable (nur im Bereich der Funktion gültig) wert = exp(x) - 20; // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(x) } // Ende der Funktion f(x) double f_strich_analyt(double x) { // Deklaration und Definition der Funktion f(x) double wert; // Lokale double-Variable (nur im Bereich der Funktion gültig) wert = exp(x); // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(x) } // Ende der Funktion f(x) int main(){ // Hauptfunktion double p[3] = {2, 2, 2}; // Deklaration des approximierten x-Wertes der Nullstelle und Start-Initialisierung const int N=10; // Anzahl der Iterationen in der Newton-Raphson Methode double h = 0.1; // Aequidistanter Abstand zwischen den x-Werten die zur numerischen Differentation benutzt werden const double p_null = log(20); // Beschreibung der ausgegebenen Groessen printf("# 0: Index i der Iteration i \n# 1: Approximierter Wert der Nullstelle p_i \n"); printf("# 2: Dreipunkte-Mittelpunkt-Formel \n# 3: Fuenfpunkte-Mittelpunkt-Formel \n"); printf("# 4: Relativer Fehler zum wirklichen Wert |(p-p_0)/p| \n# 5: Relativer Fehler zum wirklichen Wert |(p-p_1)/p| \n"); printf("# 6: Relativer Fehler zum wirklichen Wert |(p-p_2)/p| \n"); for(int i=0; i<N; ++i){ // For-Schleife der Newton-Raphson Methode printf("%3d %20.14f %20.14f %20.14f %20.14f %20.14f %20.14f \n",i, p[0], p[1], p[2], fabs((p_null-p[0])/p_null), fabs((p_null-p[1])/p_null), fabs((p_null-p[2])/p_null)); p[0] = p[0] - f(p[0])/f_strich_analyt(p[0]); // Newton-Raphson Methode mit analytischer Ableitung f'(x) p[1] = p[1] - f(p[1])/( (f(p[1]+h) - f(p[1]-h))/(2*h) ); // Newton-Raphson Methode, Dreipunkte-Mittelpunkt-Formel fuer f' p[2] = p[2] - f(p[2])/( (f(p[2]-2*h) - 8*f(p[2]-h) + 8*f(p[2]+h) - f(p[2]+2*h))/(12*h) ); // Newton-Raphson Methode, Fuenfpunkte-Mittelpunkt-Formel fuer f' } // Ende der For-Schleife der Newton-Raphson Methode printf("# Wirklicher Wert der Nullstelle: %20.14f \n", p_null); }
Das Programm erzeugt die folgende Terminalausgabe.
Man erkennt, dass das Konvergenzverhalten der 'Dreipunkte-Mittelpunkt-Formel' ein wenig schlechter als das der 'Fünfpunkte-Mittelpunkt-Formel' ist. Die berechneten Nullstellenwerte der 'Fünfpunkte-Mittelpunkt-Formel' unterscheiden sich nur unmerklich von den Werte, die sich mittels der Berechnung mit analytischer $f^\prime(x)$ Funktion ergeben.
Diese Aufgabe ist angelehnt an das Kapitel 23 "Der gedämpfte harmonische Oszillator" des Buches von Prof. Walter Greiner, Mechanik (Teil 1) [5. Auflage, 1989, siehe Seite 226- 237]. Siehe auch Vorlesungsskript von Prof. Rischke auf Seite 117- 126 http://itp.uni-frankfurt.de/~drischke/Skript_MI_WiSe2022-2023.pdf ). Wir betrachten im Folgenden den gedämpften harmonischen Oszillator am Beispiel eines reibungsfrei gelagerten Wagens (Masse=$M$) auf den eine Rückstellkraft einwirkt (die proportional zu seiner Auslenkung $x$ ist (Proportionalitätskonstante $k$)), wobei zusätzlich eine geschwindigkeitsabhängige Reibungskraft auf den Wagen einwirkt (z.B. verursacht durch den auf den Wagen einwirkende Luftwiderstand, Stokesscher Ansatz: Proportionalitätskonstante $\alpha$). Aufgrund der Rückstellkraft, besitzt das zugrundeliegende Potenzial $V(x)$ die Form einer Parabel $V(x)=\frac{k \, x^2}{2}$. Die zeitliche Entwicklung des linearen harmonischen Oszillators mit Dämpfung wird mittels der folgenden Differenzialgleichung zweiter Ordnung beschrieben (wir setzen $\omega_0^2=\frac{k}{M}$ und $\beta = \frac{\alpha}{2M}$):
$$
\begin{equation}
\ddot{x}(t) = - \omega_0^2 \, x(t) - 2 \beta \, \dot{x}(t)
\end{equation}
$$
Die Anfangsbedingungen seien zunächst noch allgemein gehalten: $x(0) = x_0 \,\, , \,\, \dot{x}(0) = v_0$. Bestimmen Sie die allgemeine Lösung der Differenzialgleichung mittels eines eigenen Jupyter Notebooks. Geben Sie dann die spezielle Lösung der Differenzialgleichung bei festgelegten Parameterwerten ($\omega_0^2=3$ und $\beta = 0.25$) und Anfangsbedingungen ($x_0 = 0$ und $v_0 = 40$) an und visualisieren Sie diese in einem x-t Diagramm. An welchem Ort befindet sich der Wagen zur Zeit $t=10$ ( $x(10)$ )? Erstellen Sie zusätzlich eine Animation, indem Sie einen der zuvor festgelegten Parameter (z.B. $x_0$, $v_0$, $\omega_0$ oder $\beta$) in einem gewissen Wertebereich verändern.
Berechnen Sie nun die allgemeine Lösung für den aperiodischen Grenzfall (kritische Dämpfung, $\beta = \omega_0$) und erstellen Sie eine weitere Animation, welche die Abbildung 2.23 im Vorlesungsskript von Prof. Rischke (Seite 124) verdeutlicht.
Das nebenstehende Jupyter Notebook A7_2.ipynb (View Notebook, Download Notebook) stellt die Musterlösung der Aufgabe dar.