Musterlösung zum Übungsblatt Nr. 6
Musterlösung zur Aufgabe 1 (6.5 Punkte)
Gegeben seien die folgenden x- und y-Werte der Stützstellen: $\vec{x} = \left( 2, 10.1, 15, 20, 40, 60, 90 \right)$ und $\vec{y} = \left( 2.1, 41, 43, 40.2, 12, 5, 17.5 \right)$. Das folgende C++ Programm berechnet das Lagrange Polynom $P_6(x)$ im Teilintervall $[a,b]=[0,100]$. Als Vorlage des Programms wurde das C++ Programm Lagrange_Polynom_1.cpp verwendet. Die wesentliche Abänderung wurde in der Kern-Gleichung der Lagrange Polynom Methode gemacht (Anweisung: "Pfp[j] += points_y[k]*Lk;") bei der nun nicht die Funktionswerte "f(points[k])", sondern die "points_y[k]"-Werte stehen.
// Musterlösung der Aufgabe 1 des Übungsblattes Nr.6 /* Lagrange Polynom von einer Liste von Stützstellen * Ausgabe zum Plotten (Gnuplot oder Python) mittels: "./a.out >> A6_1.dat" */ #include <iostream> // Ein- und Ausgabebibliothek int main(){ // Hauptfunktion double points[] = { 2, 10.1, 15, 20, 40, 60, 90 }; // Deklaration und Initialisierung der Stützstellen als double-Array double points_y[] = { 2.1, 41, 43, 40.2, 12, 5, 17.5 }; // Deklaration und Initialisierung der Stützstellen als double-Array const unsigned int N_points = sizeof(points)/sizeof(points[0]); // Anzahl der Punkte die zur Approximation verwendet werden double a = 1; // Untergrenze des x-Intervalls in dem die Ergebnisse ausgegeben werden sollen double b = 90.7; // Obergrenze des x-Intervalls in dem die Ergebnisse ausgegeben werden sollen const unsigned int N_xp = 500; // 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; // Deklaration des aktuellen x-Wertes double xp[N_xp+1]; // Deklaration der x-Ausgabe-Punkte als double-Array double Pfp[N_xp+1]; // Deklaration der Ausgabe-Punkte des approximierten Polynoms als double-Array double Lk; // Deklaration einer Variable, die fuer Zwischenrechnungen benoetigt wird printf("# x-Werte der %3d Stuetzstellen-Punkte: \n", N_points); // Beschreibung der ausgegebenen Groessen for(unsigned int k = 0; k < N_points; ++k){ // For-Schleife der Ausgabe der Stuetzstellen x-Werte printf("%10.5f",points[k]); // Ausgabe der Stuetzpunkte } // Ende for-Schleife der Ausgabe printf("\n"); // Zeilenumbruch printf("# y-Werte der %3d Stuetzstellen-Punkte: \n", N_points); // Beschreibung der ausgegebenen Groessen for(unsigned int k = 0; k < N_points; ++k){ // For-Schleife der Ausgabe der Stuetzstellen y-Werte printf("%10.5f",points_y[k]); // Ausgabe der Stuetzpunkte } // Ende for-Schleife der Ausgabe printf("\n"); // Zeilenumbruch printf("# 0: Index j \n# 1: x-Wert \n"); // Beschreibung der ausgegebenen Groessen printf("# 2: Approximierter Wert des Lagrange Polynoms P(x) \n"); // Beschreibung der ausgegebenen Groessen for(unsigned int j = 0; j <= N_xp; ++j){ // For-Schleife die ueber die einzelnen Punkte des x-Intervalls geht x = a + j*dx; // Aktueller x-Wert xp[j]=x; // Eintrag des aktuellen x-Wertes in das x-Array Pfp[j]=0; // Initialisierung des j-ten Polynom-Arrays mit 0 (Wert beim Anfang der Summenbildung) for(unsigned int k = 0; k < N_points; ++k){ // For-Schleife der Summation in der Lagrange Polynom Methode Lk=1; // Initialisierung der Produktvariable Lk mit 1 for(unsigned int i = 0; i < N_points; ++i){ // For-Schleife der Produktbildung in der Lagrange Polynom Methode if(i != k){ // Die Produktbildung soll nur fuer (i ungleich k) erfolgen Lk *= (x - points[i])/(points[k] - points[i]); // Berechnung der Lk-Werte in der Lagrange Polynom Methode } // Ende if-Bedingung } // Ende for-Schleife der Produktbildung Pfp[j] += points_y[k]*Lk; // Kern-Gleichung in der Lagrange Polynom Methode } // Ende for-Schleife der Summenbildung } // Ende der for-Schleife ueber die einzelnen Punkte des x-Intervalls for(unsigned int j = 0; j <= N_xp; ++j){ // For-Schleife der separaten Ausgabe der berechneten Werte printf("%3d %14.10f %14.10f \n",j, xp[j], Pfp[j]); // Ausgabe der berechneten Werte } // Ende for-Schleife der Ausgabe } // Ende der Hauptfunktion
Die nebenstehende Abbildung visualisiert die Ergebnisse des Programms (das zugehörige Python Skript findet man unter PythonPlot_A6_1.py).
Musterlösung zur Aufgabe 2 (7 Punkte)
Für die Vektoren $\vec{a}$ und $\vec{b}$ \[ \begin{equation} \vec{a} =\left( \begin{array}{ccc} 2.3 \\ -1.36 \\ 6.91 \end{array} \right) \quad , \quad \vec{b} =\left( \begin{array}{ccc} 10.3 \\ -4.34 \\ 5.3 \end{array} \right) \end{equation} \] wurde das Skalar- und Kreuzprodukt mittels der folgenden Ausdrücke bestimmt: \[ \begin{eqnarray} \hbox{Skalarprodukt:} && s = \vec{a} \cdot \vec{b} = \sum_{i=1}^3 a_i \, b_i \\ \hbox{Kreuzprodukt:} && \vec{c} = \vec{a} \times \vec{b} = \sum_{i,j,k=1}^3 \epsilon_{ijk} \, a_i \, b_j \, \vec{e}_k \quad , \end{eqnarray} \] wobei $\epsilon_{ijk}$ der total antisymmetrischer Tensor (auch Epsilon-Tensor bzw. Levi-Civita-Symbol) ist. Das folgende C++ Programm stellt den Quelltext der Musterlösung dar.
/* Musterlösung der Aufgabe 2 des Übungsblattes Nr.6 * Skalarprodukt und Kreuzprodukt zweier Vektoren * deklariert mittels zweier eindimensionaler Arrays */ #include <iostream> // Ein- und Ausgabebibliothek int main(){ // Hauptfunktion double a[] = { 2.3, -1.36, 6.9}; // Definition des Vektors a mittels eines double-Arrays double b[] = {10.3, -4.34, 5.3}; // Definition des Vektors b mittels eines double-Arrays double c_1[3] = {0, 0, 0}; // Deklaration des Vektors c mittels eines double-Arrays (c = a x b) double c_2[3] = {0, 0, 0}; // Deklaration des Vektors c mittels eines double-Arrays (c = a x b) int epsilon[3][3][3]; // Total antisymmetrischer Tensor (Levi-Civita-Symbol, Epsilon-Tensor) double skal_prod_1=0; // Deklaration der Variable fuer das Skalarprodukt a*b mittels Index-Zugriff double skal_prod_2=0; // Deklaration der Variable fuer das Skalarprodukt a*b mittels Zeiger-Zugriff double* pa = a; // Definition des Zeigers auf das double-Array a double* pb = b; // Definition des Zeigers auf das double-Array b double* pc_2 = c_2; // Definition des Zeigers auf das double-Array c_2 // Skalarprodukt a*b for(int i=0; i<3; ++i){ skal_prod_1 = skal_prod_1 + a[i]*b[i]; skal_prod_2 = skal_prod_2 + *(pa+i) * *(pb+i); } printf("Das Skalarprodukt mittels Index-Zugriff: a*b = %8.5f \n", skal_prod_1); printf("Das Skalarprodukt mittels Zeiger-Zugriff: a*b = %8.5f \n", skal_prod_2); // Initialisierung der Werte für den Epsilon-Tensor for(int i=0; i<3; ++i){ for(int j=0; j<3; ++j){ for(int k=0; k<3; ++k){ if (i==j || i==k || j==k){ epsilon[i][j][k] = 0; } else if ( (i==0 && j==1 && k==2) || (i==1 && j==2 && k==0) || (i==2 && j==0 && k==1)){ epsilon[i][j][k] = 1; } else { epsilon[i][j][k] = -1; } } } } // Kreuzprodukt c = a x b for(int i=0; i<3; ++i){ for(int j=0; j<3; ++j){ for(int k=0; k<3; ++k){ c_1[k] = c_1[k] + epsilon[i][j][k] * a[i] * b[j]; *(pc_2+k) = *(pc_2+k) + epsilon[i][j][k] * *(pa+i) * *(pb+j); } } } printf("Das Kreuzprodukt mittels Index-Zugriff: a x b = ( %8.5f , %8.5f, %8.5f ) \n", c_1[0], c_1[1], c_1[2]); printf("Das Kreuzprodukt mittels Zeiger-Zugriff: a x b = ( %8.5f , %8.5f, %8.5f ) \n", c_2[0], c_2[1], c_2[2]); }
Es wurden zunächst die angegebenen Vektoren $\vec{a}$ und $\vec{b}$ als eindimensionale C++ Arrays und auch zwei entsprechende Zeigervariablen (double* pa = a; und double* pb = b; ) definiert. Dann wurde das Skalarprodukt mittels eines Indexzugriffs ( skal_prod_1 = skal_prod_1 + a[i]*b[i]; ) bzw. Zeigerzugriffs ( skal_prod_2 = skal_prod_2 + *(pa+i) * *(pb+i); ) realisiert. Der schon weiter oben deklarierte Epsilon-Tensor $\epsilon_{ijk}$ wurde dann mittels dreier geschachtelter for-Schleifen initialisiert. Die darauf folgende Berechnung des Kreuzproduktes wurde ebenfalls mittels eines Indexzugriffs ( c_1[k] = c_1[k] + epsilon[i][j][k] * a[i] * b[j]; ) bzw. Zeigerzugriffs ( *(pc_2+k) = *(pc_2+k) + epsilon[i][j][k] * *(pa+i) * *(pb+j); ) realisiert. Die nebenstehende Abbildung zeigt die Terminalausgabe des Programms.
Musterlösung zur Aufgabe 3 (6.5 Punkte)
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 3 des Übungsblattes Nr.6 #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.