Numerische Differentiation
In diesem Unterpunkt werden wir die im Jupyter Notebook Das Computeralgebrasystem: Python Jupyter + SymPy ('Jupyter_sympy.ipynb' , View Notebook, Download Notebook) hergeleiteten Differentiationsregeln der numerischen Mathematik in einem C++ Programm anwenden. Im Speziellen werden wir die numerische Ableitung der Funktion $f(x)= 10 \, e^{-x/10} \cdot \hbox{sin}(x)$ mittels der unterschiedlichen Differentiationsformel in einem C++ Programm bestimmen und die Ergebnisse mit der wirklichen, analytisch bestimmbaren Ableitung $f^\prime(x)$ vergleichen.
Die mathematische Definition der Ableitung einer Funktion $f(x)$ an der Stelle $x_0$, der sogenannte Differentialquotient, benutzt die folgende Grenzwert-Formulierung ($h \to 0$):
\[
\begin{equation}
f^\prime(x_0) = \lim \limits_{h \to 0} \frac{f(x_0+h) - f(x_0)}{h} \quad.
\end{equation}
\]
In der numerischen Mathematik hingegen werden mehrere Differentiationsregeln formuliert, die diesen wirklichen Wert der Ableitung approximieren. Die analytische Herleitung dieser Regeln benutzt die Methode der Lagrange Polynome (siehe Anwendungsbeispiel: Interpolation und Polynomapproximation).
Die Differentiationsregeln der numerischen Mathematik
Die im vorigen Unterpunkt hergeleiteten Differentiationsregeln der numerischen Mathematik (siehe Das Computeralgebrasystem: Python Jupyter + SymPy) sind in der folgenden Box nochmals zusammenfassend dargestellt:
\[
\begin{eqnarray}
\hbox{Zweipunkteformel (n=1):}\,\, & \,\,
f^\prime(x_0) \approx \frac{1}{h} \left[ f(x_0 + h) - f(x_0) \right] & \, , \, \hbox{Stützstellen: }\, (x_0, x_0 + h) \\
\hbox{Dreipunkte-Endpunkt-Formel (n=2):}\,\, & \,\,
f^\prime(x_0) \approx \frac{1}{2h} \left[ - 3\,f(x_0) + 4 \, f(x_0 + h) - f(x_0 + 2 \, h) \right] & \, , \, \hbox{Stützstellen: }\, (x_0, x_0 + h, x_0 + 2h) \\
\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{Stützstellen: }\, (x_0 - h, x_0, x_0 + h) \\
\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] & \, , \, \hbox{Stützstellen: }\, (x_0 - 2h, x_0 - h, x_0, x_0 + h, x_0 + 2h)
\end{eqnarray}
\]
In dem folgenden C++ Programm wird die Ableitung der Funktion $f(x)= 10 \, e^{-x/10} \cdot \hbox{sin}(x)$ an der Stelle $x_0 = 3$ numerisch approximativ berechnet.
Numerical_Differentiation_1.cpp
/* Berechnung der Ableitung f'(x) einer Funktion f(x) * mittels unterschiedlich genauer Approximationen * (Zweipunkteformel, Dreipunkte-Endpunkt-Formel, Dreipunkte-Mittelpunkt-Formel und Fünfpunkte-Mittelpunkt-Formel) */ #include <stdio.h> // Standard Input- und Output Bibliothek in C, z.B. printf(...) #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) double f(double x){ // Deklaration und Definition der Funktion f(x) double wert; wert=10*sin(x)*exp(-x/10); // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(x) } // Ende der Funktion f(x) double f_strich(double x){ // Deklaration und Definition der Ableitung f'(x) der Funktion f(x) double wert; // (f'(x) wird nur zum Vergleich mit den Ergebnissen benoetigt) wert=10*cos(x)*exp(-x/10) - sin(x)*exp(-x/10); // Eigentliche Definition der Ableitung return wert; // Rueckgabewert der Funktion f'(x) } // Ende der Funktion f'(x) int main(){ // Hauptfunktion double x = 3; // Aktueller x-Wert an dem die Ableitung berechnet wird double h=0.01; // Aequidistanter Abstand zwischen den x-Werten die zur numerischen Differentation benutzt werden double f1_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Zweipunkteformel) double f3a_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Dreipunkte-Endpunkt-Formel) double f3b_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Dreipunkte-Mittelpunkt-Formel) double f5_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Fünfpunkte-Mittelpunkt-Formel) f1_strich_p = (f(x+h) - f(x))/h; // Zweipunkteformel f3a_strich_p = (-3*f(x) + 4*f(x+h) - f(x+2*h))/(2*h); // Dreipunkte-Endpunkt-Formel f3b_strich_p = (f(x+h) - f(x-h))/(2*h); // Dreipunkte-Mittelpunkt-Formel f5_strich_p = (f(x-2*h) - 8*f(x-h) + 8*f(x+h) - f(x+2*h))/(12*h); // Fünfpunkte-Mittelpunkt-Formel printf("h-Wert: h = %14.10f \n", h); printf("Wirklicher Wert: f'(x=%5.3f) = %14.10f \n", x, f_strich(x)); // Ausgabe der berechneten Werte printf("Zweipunkteformel: f'(x=%5.3f) = %14.10f , absoluter Fehler: %14.10f \n", x, f1_strich_p, f_strich(x) - f1_strich_p); // Ausgabe der berechneten Werte printf("Dreipunkte-Endpunkt-Formel: f'(x=%5.3f) = %14.10f , absoluter Fehler: %14.10f \n", x, f3a_strich_p, f_strich(x) - f3a_strich_p); // Ausgabe der berechneten Werte printf("Dreipunkte-Mittelpunkt-Formel: f'(x=%5.3f) = %14.10f , absoluter Fehler: %14.10f \n", x, f3b_strich_p, f_strich(x) - f3b_strich_p); // Ausgabe der berechneten Werte printf("Fuenfpunkte-Mittelpunkt-Formel: f'(x=%5.3f) = %14.10f , absoluter Fehler: %14.10f \n", x, f5_strich_p, f_strich(x) - f5_strich_p); // Ausgabe der berechneten Werte } // Ende der Hauptfunktion
Es wird hierbei ein fester Wert von $h=0.01$ angenommen und der Wert der Ableitung $f^\prime(x_0)$ mittels der Zweipunkteformel, Dreipunkte-Endpunkt-Formel, Dreipunkte-Mittelpunkt-Formel und Fünfpunkte-Mittelpunkt-Formel approximiert. Die numerischen Ergebnisse werden dann mit dem analytischen Wert der Ableitung verglichen. Die nebenstehende Abbildung zeigt die Terminalausgabe des Programms.
Nun wollen wir uns die Ableitung der Funktion in einem Teilintervall $[a,b]$ ausgeben lassen und die numerischen Lösungen des C++ Programms für unterschiedliche Werte von $h$ mittels eines Python Skriptes visualisieren. Das folgende C++ Programm berechnet $300$ Werte von $f^\prime(x)$ im Teilintervall $[0.5, 10]$ und benutzt dabei drei unterschiedliche $h$-Werte ($h=0.5, 0.1, 0.01$).
Numerical_Differentiation_2.cpp
/* Berechnung der Ableitung f'(x) einer Funktion f(x) * mittels unterschiedlich genauer Approximationen * (Zweipunkteformel, Dreipunkte-Endpunkt-Formel, Dreipunkte-Mittelpunkt-Formel und Fünfpunkte-Mittelpunkt-Formel) * Ausgabe fuer unterschiedliche x-Werte in [a,b] und 3 unterschiedliche h-Werte * Ausgabe zum Plotten mittels Python: "./a.out > Numerical_Differentiation_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(double x){ // Deklaration und Definition der Funktion f(x) double wert; wert = 10*sin(x)*exp(-x/10); // Eigentliche Definition der Funktion return wert; // Rueckgabewert der Funktion f(x) } // Ende der Funktion f(x) double f_strich(double x){ // Deklaration und Definition der Ableitung f'(x) der Funktion f(x) double wert; // (f'(x) wird nur zum Vergleich mit den Ergebnissen benoetigt) wert = 10*cos(x)*exp(-x/10) - sin(x)*exp(-x/10); // Eigentliche Definition der Ableitung return wert; // Rueckgabewert der Funktion f'(x) } // Ende der Funktion f'(x) int main(){ // Hauptfunktion double a = 0.5; // Untergrenze des Intervalls [a,b] in dem f(x)=0 gelten soll double b = 10; // Obergrenze des Intervalls [a,b] in dem f(x)=0 gelten soll int N_xp = 300; // Anzahl der Punkte in die das x-Intervall aufgeteilt wird double dx = (b - a)/N_xp; // Abstand dx zwischen den aequidistanten Punkten des x-Intervalls double x = a - dx; // Aktueller x-Wert double h_list[3]={0.5, 0.1, 0.01}; // Aequidistanter Abstand zwischen den x-Werten die zur numerischen Differentation benutzt werden double h; // Aktueller h-Wert double f1_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Zweipunkteformel) double f3a_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Dreipunkte-Endpunkt-Formel) double f3b_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (3Dreipunkte-Mittelpunkt-Formel) double f5_strich_p; // Deklaration der f'(x)-Ausgabe-Punkte (Fünfpunkte-Mittelpunkt-Formel) printf("# Benutzte h-Werte der Simulation: \n"); // Beschreibung der ausgegebenen Groessen for(int k = 0; k < 3; ++k){ // For-Schleife der Ausgabe der Stuetzstellen x-Werte printf("%10.5f",h_list[k]); // Ausgabe der unterschiedlichen h-Werte } // Ende for-Schleife der Ausgabe printf("\n"); printf("# 0: Index i \n# 1: x-Wert \n# 2: Wirklicher Wert von f'(x) \n"); // Beschreibung der ausgegebenen Groessen printf("# h = %8.7f : \n", h_list[0]); // 0-ter h-Wert printf("# 3: Zweipunkteformel f'(x) \n# 4: Dreipunkte-Endpunkt-Formel \n"); // Beschreibung der ausgegebenen Groessen printf("# 5: Dreipunkte-Mittelpunkt-Formel \n# 6: Fünfpunkte-Mittelpunkt-Formel \n"); // Beschreibung der ausgegebenen Groessen printf("# h = %8.7f : \n", h_list[1]); // 1-ter h-Wert printf("# 7: Zweipunkteformel f'(x) \n# 8: Dreipunkte-Endpunkt-Formel \n"); // Beschreibung der ausgegebenen Groessen printf("# 9: Dreipunkte-Mittelpunkt-Formel \n# 10: Fünfpunkte-Mittelpunkt-Formel \n"); // Beschreibung der ausgegebenen Groessen printf("# h = %8.7f : \n", h_list[2]); // 2-ter h-Wert printf("# 11: Zweipunkteformel f'(x) \n# 12: Dreipunkte-Endpunkt-Formel \n"); // Beschreibung der ausgegebenen Groessen printf("# 13: Dreipunkte-Mittelpunkt-Formel \n# 14: Fünfpunkte-Mittelpunkt-Formel \n"); // Beschreibung der ausgegebenen Groessen for(int i =0; i < N_xp; ++i){ // for-Schleife die ueber die einzelnen Punkte des x-Intervalls geht x=x+dx; // Aktueller x-Wert printf("%3d %14.10f %14.10f",i, x, f_strich(x)); // Ausgaben for(int k = 0; k < 3; ++k){ // For-Schleife die ueber die unterschiedlichen h-Werte h = h_list[k]; // Uebergabe des aktuellen h-Wertes f1_strich_p = (f(x+h) - f(x))/h; // Zweipunkteformel f3a_strich_p = (-3*f(x) + 4*f(x+h) - f(x+2*h))/(2*h); // Dreipunkte-Endpunkt-Formel f3b_strich_p = (f(x+h) - f(x-h))/(2*h); // Dreipunkte-Mittelpunkt-Formel f5_strich_p = (f(x-2*h) - 8*f(x-h) + 8*f(x+h) - f(x+2*h))/(12*h); // Fünfpunkte-Mittelpunkt-Formel printf(" %14.10f %14.10f %14.10f %14.10f", f1_strich_p, f3a_strich_p, f3b_strich_p, f5_strich_p); // Ausgaben der berechneten Groessen } // Ende For-Schleife ueber die unterschiedlichen h-Werte printf("\n"); // Zeilenumbruch } // Ende for-Schleife ueber die einzelnen Punkte des x-Intervalls } // Ende der Hauptfunktion
Das Programm gibt die folgende Tabelle im Terminal aus:
Wir leiten die berechneten Daten wieder in eine Datei "Numerical_Differentiation_2.dat" um und führen das folgende Python Skript aus:
PythonPlot_Numerical_Differentiation_2.py
# Python Programm zum Plotten der berechneten Daten von Numerical_Differentiation_2.cpp import matplotlib # Python Bibliothek zum Plotten (siehe https://matplotlib.org/ ) 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/ ) import matplotlib.gridspec as gridspec # Mehrere Bilder in einem Plot anordnen # Bildabmessungen usw. params = { 'figure.figsize' : [14,10], 'text.usetex' : True, 'axes.titlesize' : 14, 'axes.labelsize' : 16, 'xtick.labelsize' : 14 , 'ytick.labelsize' : 14 } matplotlib.rcParams.update(params) h_list = np.genfromtxt("./Numerical_Differentiation_2.dat", max_rows=1) # Einlesen der benutzten x-Werte der Stuetzstellen data = np.genfromtxt("./Numerical_Differentiation_2.dat", skip_header=20) # Einlesen der berechneten Daten von Numerical_Differentiation_2.cpp fig = plt.figure() # Hauptbild gs = gridspec.GridSpec(2, 2, width_ratios=[1,1], wspace=0.3, hspace=0.3) # Anordnung der vier Unterbilder ax1 = plt.subplot(gs[0]) ax2 = plt.subplot(gs[1]) ax3 = plt.subplot(gs[2]) ax4 = plt.subplot(gs[3]) ax1.set_title(r'Zweipunkteformel') # Titel der Abbildung in ax1 ax1.set_xlabel(r"$\rm x$") ax1.set_ylabel(r"$\rm f'(x)$") ax2.set_title(r'Dreipunkte-Endpunkt-Formel') # Titel der Abbildung in ax2 ax2.set_xlabel(r"$\rm x$") ax2.set_ylabel(r"$\rm f'(x)$") ax3.set_title(r'Dreipunkte-Mittelpunkt-Formel') # Titel der Abbildung in ax3 ax3.set_xlabel(r"$\rm x$") ax3.set_ylabel(r"$\rm f'(x)$") ax4.set_title(r'Fünfpunkte-Mittelpunkt-Formel') # Titel der Abbildung in ax4 ax4.set_xlabel(r"$\rm x$") ax4.set_ylabel(r"$\rm f'(x)$") l_width=0.8 # Festlegung der Plot-Liniendicke alp=0.7 # Festlegung der Transparenz der Kurven label_0=r'$\rm h='+'{0:.3f}'.format(h_list[0])+'$' # Plot-Labels fuer die unterschiedlichen h-Werte label_1=r'$\rm h='+'{0:.3f}'.format(h_list[1])+'$' label_2=r'$\rm h='+'{0:.3f}'.format(h_list[2])+'$' ax1.plot(data[:,1],data[:,2], color="black", linewidth=l_width, linestyle='-.', alpha=alp, label=r'$\rm Wirklicher \, Wert$') # Plotten auf ax1 ax1.plot(data[:,1],data[:,3], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=label_0) ax1.plot(data[:,1],data[:,7], color="red", linewidth=l_width, linestyle='-', alpha=alp, label=label_1) ax1.plot(data[:,1],data[:,11], color="green", linewidth=l_width, linestyle='-', alpha=alp, label=label_2) ax2.plot(data[:,1],data[:,2], color="black", linewidth=l_width, linestyle=':', alpha=alp, label=r'$\rm Wirklicher \, Wert$') # Plotten auf ax2 ax2.plot(data[:,1],data[:,4], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=label_0) ax2.plot(data[:,1],data[:,8], color="red", linewidth=l_width, linestyle='-', alpha=alp, label=label_1) ax2.plot(data[:,1],data[:,12], color="green", linewidth=l_width, linestyle='-', alpha=alp, label=label_2) ax3.plot(data[:,1],data[:,2], color="black", linewidth=l_width, linestyle=':', alpha=alp, label=r'$\rm Wirklicher \, Wert$') # Plotten auf ax3 ax3.plot(data[:,1],data[:,5], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=label_0) ax3.plot(data[:,1],data[:,8], color="red", linewidth=l_width, linestyle='-', alpha=alp, label=label_1) ax3.plot(data[:,1],data[:,12], color="green", linewidth=l_width, linestyle='-', alpha=alp, label=label_2) ax4.plot(data[:,1],data[:,2], color="black", linewidth=l_width, linestyle=':', alpha=alp, label=r'$\rm Wirklicher \, Wert$') # Plotten auf ax4 ax4.plot(data[:,1],data[:,6], color="blue", linewidth=l_width, linestyle='-', alpha=alp, label=label_0) ax4.plot(data[:,1],data[:,9], color="red", linewidth=l_width, linestyle='-', alpha=alp, label=label_1) ax4.plot(data[:,1],data[:,13], color="green", linewidth=l_width, linestyle='-', alpha=alp, label=label_2) ax1.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax1 ax2.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax2 ax3.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax3 ax4.legend(frameon=True, loc="lower right",fontsize=10) # Anordnung der Legende auf ax4 plt.savefig("Numerical_Differentiation_2.png", dpi=400, bbox_inches="tight", pad_inches=0.05, format="png") # Speichern der Abbildung als Bild plt.show() # Zusaetzliches Darstellen der Abbildung in einem separaten Fenster
Das Python Skript produziert das folgende Bild:
