diff --git a/sose2020/num/uebungen/num10.pdf b/sose2020/num/uebungen/num10.pdf new file mode 100644 index 0000000..cf4e0c7 Binary files /dev/null and b/sose2020/num/uebungen/num10.pdf differ diff --git a/sose2020/num/uebungen/num10.tex b/sose2020/num/uebungen/num10.tex new file mode 100644 index 0000000..7cfe58b --- /dev/null +++ b/sose2020/num/uebungen/num10.tex @@ -0,0 +1,262 @@ +\documentclass[uebung]{../../../lecture} + +\title{Einführung in die Numerik: Übungsblatt 10} +\author{Leon Burgard, Christian Merten} + +\begin{document} + +\punkte + +\begin{aufgabe} + Beh.: Seien $n$ paarweise verschiedene Stützstellen $\{x_0, \ldots, x_{n-1}\} $ gegeben + und eine Permutation derselben $\{\tilde{x}_0, \ldots, \tilde{x}_{n-1}\} $. Dann gilt + \[ + f[x_0, \ldots, x_{n-1}] = f[\tilde{x}_0, \ldots, \tilde{x}_{n-1}] + .\] + \begin{proof} + Es gilt nach VL mit der Newtondarstellung für das Interpolationspolynom + zu den Stützstellen $x_0, \ldots, x_{n-1}$: + \begin{salign*} + p(x) &= \sum_{i=0}^{n-1} y[x_0, \ldots, x_i] N_i(x) \\ + &= \sum_{i=0}^{n-2} y[x_0, \ldots, x_i] N_i(x) + y[x_0, \ldots, x_{n-1}] N_{n-1}(x) \\ + &= \mathcal{O}(x^{n-2}) + y[x_0, \ldots, x_{n-1}] \prod_{i=0}^{n-2} (x - x_i) \\ + &= \mathcal{O}(x^{n-2}) + y[x_0, \ldots, x_{n-1}] \left(x^{n-1} + \mathcal{O}(x^{n-2})\right) \\ + &= y[x_0, \ldots, x_{n-1}] x^{n-1} + \mathcal{O}(x^{n-2}) + .\end{salign*} + Der Leitkoeffizient des Interpolationspolynoms in der Monombasis ist also + $y[x_0, \ldots, x_{n-1}]$. Dieser ist unabhängig von der Reihenfolge der Stützstellen. Damit + folgt die Behauptung. + \end{proof} +\end{aufgabe} + +\begin{aufgabe} + Beh.: Für $N \ge 2 \pi 10^{3}$ gilt $\displaystyle \max_{0 \le x \le 1} |f(x) - s(x)| < 10^{-12}$. + \begin{proof} + Es ist $f \in C^{4}([0, 1])$. Dann gilt nach VL + \[ + \delta \coloneqq \max_{0 \le x \le 1} |f(x) - s(x)| \le h^{4} \max_{0 \le x \le 1} |f^{(4)}(x)| + .\] Mit $f(x) = \sin(2\pi x)$ folgt sofort + \[ + f^{(4)}(x) = 16 \pi^{4} \sin(2 \pi x) \quad \text{also}\quad + \max_{0 \le x \le 1} |f^{(4)}(x)| = 16 \pi^{4} + .\] Mit $h = \frac{1}{N}$ ergibt sich + \[ + \delta \le 16 \frac{\pi^{4}}{N^{4}} \implies N \ge 2 \pi \sqrt[4]{\delta } = 2 \pi 10^{3} + .\] + \end{proof} +\end{aufgabe} + +\begin{aufgabe} + Beh.: Das komplexe trigonometrische Interpolationspolynom ist gegeben als + \[ + t ^{*}(x) = \frac{1}{2} - \frac{1}{4} e^{ix} - \frac{1}{4} e^{3ix} + .\] + \begin{proof} + Die Stützstellen sind als $x_j = \frac{2 \pi j}{4}$, $j = 0, \ldots, 3$ gegeben. Damit folgt + \begin{salign*} + f(x_0) &= f(0) = \min \{0, 2\} = 0 \\ + f(x_1) &= f\left(\frac{\pi}{2}\right) = \min \left\{ \frac{1}{2}, \frac{3}{2} \right\} = \frac{1}{2} \\ + f(x_2) &= f(\pi) = \min \{1, 1\} = 1 \\ + f(x_3) &= f\left( \frac{3}{2}\pi \right) = \min \left\{ \frac{3}{2}, \frac{1}{2} \right\} = \frac{1}{2} + .\end{salign*} + Die Interpolationsbedingung ist erfüllt, denn + \begin{salign*} + t ^{*}(x_0) &= \frac{1}{2} - \frac{1}{4} - \frac{1}{4} = 0 = f(x_0) \\ + t ^{*}(x_1) &= \frac{1}{2} - \frac{1}{4} \underbrace{e^{i \frac{\pi}{2}}}_{= i} - \frac{1}{4} + \underbrace{e^{i \frac{3}{2} \pi}}_{= -i} = \frac{1}{2} = f(x_1) \\ + t ^{*}(x_2) &= \frac{1}{2} - \frac{1}{4} \underbrace{e^{i \pi}}_{= -1} - \frac{1}{4} + \underbrace{e^{i \pi}}_{= -1} = 1 = f(x_2) \\ + t ^{*}(x_3) &= \frac{1}{2} - \frac{1}{4} \underbrace{e^{i \frac{3}{2} \pi}}_{= i} - \frac{1}{4} + \underbrace{e^{i \frac{3}{2} \pi}}_{= -i} = \frac{1}{2} = f(x_3) + .\end{salign*} + Aus der Eindeutigkeit des komplexen trigonometrischen Interpolationspolynoms folgt die Behauptung. + \end{proof} +\end{aufgabe} + +\begin{aufgabe} + \begin{enumerate}[a)] + \item + \begin{enumerate}[(1)] + \item Es ist + \[ + f_1''(x) = 6x - 14 \implies f_1''(0) = -14 \neq 0 + .\] Also erfüllt $f_1$ nicht die natürlichen Randbedingungen, also $f_1 \not\in S(X)$. + \item Es gilt für $0 \le x < 1$: + \begin{salign*} + f_2(x) &= -\frac{1}{2} x^{3} \xrightarrow{x \to 1} -\frac{1}{2}\\ + f_2'(x) &= -\frac{3}{2} x^2 \xrightarrow{x \to 1} - \frac{3}{2} \\ + f_2''(x) &= - 3x \xrightarrow{x \to 1} -3 \text{ und } f_2''(0) = 0 + \intertext{Für $1 \le x \le 2$ gilt:} + f_2(x) &= (x-1)^{3} -\frac{1}{2} x^{3} \implies f_2(1) = -\frac{1}{2} \\ + f_2'(x) &= 3(x-1)^2 - \frac{3}{2} x^2 \implies f_2'(1) = -\frac{3}{2} \\ + f_2''(x) &= 3x - 6 \implies f_2''(1) = -3 \text{ und } f_2''(2) = 0 + .\end{salign*} + $f_2$ ist auf beiden Teilintervallen ein Polynom von Grad $3$ und damit + auf den Teilintervallen beliebig oft stetig differenzierbar. Außerdem + ist $f_2$ $2$ mal stetig differenzierbar an der Stelle $1$, also insgesamt + $f_2 \in C^{2}([0, 2])$. Die natürlichen Randbedingungen sind außerdem erfüllt, also + folgt $f_2 \in S(X)$. + \item Es ist + \[ + f_3''(x) = 6x - 2 \implies f_3''(0) = -2 \neq 0 + .\] Also erfüllt $f_3$ nicht die natürlichen Randbedingungen, also $f_3 \not\in S(X)$. + \end{enumerate} + \item Der interpolierende Spline $s$ von $f(x) = x^{3}$ folgt mit der Darstellung der VL direkt + als + \[ + s(x) = \begin{cases} + 1 + 4x + 4,5 x^2 + 1,5 x^{3} & x \in [0, 1) \\ + 8 + 8,5(x-1) - 1,5(x-1)^{3} & x \in [1,2] + \end{cases} + .\] + \end{enumerate} +\end{aufgabe} + +\begin{aufgabe} + Auszüge aus \textit{splines.cc}: + \begin{lstlisting}[language=C++, title=Schneller Löser für tridiagonale Matrizen, captionpos=b] +template +void solveTriDiag(hdnum::DenseMatrix &A, std::vector &x, std::vector &b) { + int N = b.size(); + x[0] = A[0][0]; + // LU Zerlegung + REAL l; + for (int j=1; j=0; j--) { + x[j] = (b[j] - A[j][j+1]*x[j+1])/x[j]; + } +}\end{lstlisting} +Implementation in einer Klasse \lstinline{CubicSpline}. Die Funktion \lstinline{getCubicSpline} +entspricht dem ersten Konstruktor. +\begin{lstlisting}[language=C++, title=Konstruktion und Auswertung eines kubischen Splines, captionpos=b] +template +class CubicSpline { + public: + // erstelle einen kubischen spline mit vorgegebenen stuetzstellen + // und werten + CubicSpline(std::vector xs, std::vector ys) { + calculateCoefficients(xs, ys); + } + + // erstelle einen kubischen spline mit vorgegebenen stuetzstellen + // und einer zu interpolierenden funktion + CubicSpline(std::vector xs, REAL(*f)(REAL)) { + int N = xs.size(); + std::vector ys(N); + for (int i=0; i xs(N+1); + std::vector ys(N+1); + for (int i=0; i<=N; i++) { + xs[i] = a + (1.0*i)/N*(b-a); + ys[i] = f(xs[i]); + } + calculateCoefficients(xs, ys); + } + + // werte kubischen spline an vorgegebener stelle aus + REAL evaluate(REAL x) { + for (int i=1; i x_s[i] && i < x_s.size() - 1) { + continue; + } else { + return a_0[i-1] + a_1[i-1] *(x - x_s[i]) + a_2[i-1]*std::pow(x - x_s[i], 2) + a_3[i-1]*std::pow(x-x_s[i], 3); + } + } + return 0; + } + + // gebe alle interpolations polynome aus + void print() { + for(int i=0; i x_s; + // koeffizienten + std::vector a_0; + std::vector a_1; + std::vector a_2; + std::vector a_3; + + void calculateCoefficients(std::vector xs, std::vector ys) { + // stuetzstellen from x_0 ... to x_n + int n = xs.size()-1; + // copy stuetzstellen + x_s = std::vector(n+1); + x_s = xs; + // setup (n-1)x(n-1) matrix for a_2 + hdnum::DenseMatrix A(n-1,n-1); + std::vector b(n-1); + std::vector x(n-1); + REAL h; // h_i + REAL h1; // h_{i+1} + // setup LGS for a_2 + // A has tridiagonal structure + for (int i = 1; i 1) { + A[i-1][i-2] = h; + } if (i < n-1) { + A[i-1][i] = h1; + } + A[i-1][i-1] = 2 * (h + h1); + } + // initialize vectors + a_0 = std::vector(n); + a_1 = std::vector(n); + a_2 = std::vector(n); + a_3 = std::vector(n); + solveTriDiag(A,a_2,b); + // natuerliche randbedingung + a_2[n-1] = 0; + // berechne restliche koeffizienten + for (int i = 1; i<=n; i++) { + h = xs[i] - xs[i-1]; // h_i + h1 = xs[i+1] - xs[i]; // h_{i+1} + a_0[i-1] = ys[i]; + if (i == 1) { // a_2[-1] = 0 + a_1[i-1] = (ys[i] - ys[i-1])/h + (h/3)*(2*a_2[i-1]); + a_3[i-1] = (a_2[i-1])/(3*h); + } else { + a_1[i-1] = (ys[i] - ys[i-1])/h + h/3*(2*a_2[i-1] + a_2[i-2]); + a_3[i-1] = (a_2[i-1] - a_2[i-2])/(3 * h); + } + } + } +};\end{lstlisting} +\begin{figure}[h] + \centering + \begin{tikzpicture} + \begin{axis}[xtick=\empty, ytick=\empty] + \addplot[purple] table {saurier.dat}; + \end{axis} + \end{tikzpicture} + \caption{Rekonstruktion des Sauriers} + \label{fig:} +\end{figure} + +\end{aufgabe} + +\end{document} diff --git a/sose2020/num/uebungen/splines.cc b/sose2020/num/uebungen/splines.cc new file mode 100644 index 0000000..c65d077 --- /dev/null +++ b/sose2020/num/uebungen/splines.cc @@ -0,0 +1,161 @@ +#include +#include "hdnum.hh" // hdnum header + +template +void solveTriDiag(hdnum::DenseMatrix &A, std::vector &x, std::vector &b) { + int N = b.size(); + x[0] = A[0][0]; + // LU Zerlegung + REAL l; + for (int j=1; j=0; j--) { + x[j] = (b[j] - A[j][j+1]*x[j+1])/x[j]; + } +} + +template +class CubicSpline { + public: + // erstelle einen kubischen spline mit vorgegebenen stuetzstellen + // und werten + CubicSpline(std::vector xs, std::vector ys) { + calculateCoefficients(xs, ys); + } + + // erstelle einen kubischen spline mit vorgegebenen stuetzstellen + // und einer zu interpolierenden funktion + CubicSpline(std::vector xs, REAL(*f)(REAL)) { + int N = xs.size(); + std::vector ys(N); + for (int i=0; i xs(N+1); + std::vector ys(N+1); + for (int i=0; i<=N; i++) { + xs[i] = a + (1.0*i)/N*(b-a); + ys[i] = f(xs[i]); + } + calculateCoefficients(xs, ys); + } + + // werte kubischen spline an vorgegebener stelle aus + REAL evaluate(REAL x) { + for (int i=1; i x_s[i] && i < x_s.size() - 1) { + continue; + } else { + return a_0[i-1] + a_1[i-1] *(x - x_s[i]) + a_2[i-1]*std::pow(x - x_s[i], 2) + a_3[i-1]*std::pow(x-x_s[i], 3); + } + } + return 0; + } + + // gebe alle interpolations polynome aus + void print() { + for(int i=0; i x_s; + // koeffizienten + std::vector a_0; + std::vector a_1; + std::vector a_2; + std::vector a_3; + + void calculateCoefficients(std::vector xs, std::vector ys) { + // stuetzstellen from x_0 ... to x_n + int n = xs.size()-1; + // copy stuetzstellen + x_s = std::vector(n+1); + x_s = xs; + // setup (n-1)x(n-1) matrix for a_2 + hdnum::DenseMatrix A(n-1,n-1); + std::vector b(n-1); + std::vector x(n-1); + REAL h; // h_i + REAL h1; // h_{i+1} + // setup LGS for a_2 + // A has tridiagonal structure + for (int i = 1; i 1) { + A[i-1][i-2] = h; + } if (i < n-1) { + A[i-1][i] = h1; + } + A[i-1][i-1] = 2 * (h + h1); + } + // initialize vectors + a_0 = std::vector(n); + a_1 = std::vector(n); + a_2 = std::vector(n); + a_3 = std::vector(n); + solveTriDiag(A,a_2,b); + // natuerliche randbedingung + a_2[n-1] = 0; + // berechne restliche koeffizienten + for (int i = 1; i<=n; i++) { + h = xs[i] - xs[i-1]; // h_i + h1 = xs[i+1] - xs[i]; // h_{i+1} + a_0[i-1] = ys[i]; + if (i == 1) { // a_2[-1] = 0 + a_1[i-1] = (ys[i] - ys[i-1])/h + (h/3)*(2*a_2[i-1]); + a_3[i-1] = (a_2[i-1])/(3*h); + } else { + a_1[i-1] = (ys[i] - ys[i-1])/h + h/3*(2*a_2[i-1] + a_2[i-2]); + a_3[i-1] = (a_2[i-1] - a_2[i-2])/(3 * h); + } + } + } +}; + +double f1(double x) { + return std::pow(x,3); +} + +int main() { + // std::vector xs = {0, 1, 2}; + // CubicSpline phi(xs, f1); + // phi.print(); + + // saurier fundstellen + std::vector xs = {3.75, 3.75, 2.25, 1.25, 0.25, 0.75, 4.00, 5.50, 6.25, 9.00, 12.5, 12.75, 12.25}; + std::vector ys = {0.25, 1.50, 3.25, 4.25, 4.65, 4.83, 2.50, 3.25, 3.65, 3.00, 1.25, 2.150, 3.500}; + + std::vector ts(13); + ts[0] = 0; + for (int i = 1; i<=12; i++) { + ts[i] = ts[i-1] + std::sqrt(std::pow(xs[i] - xs[i-1], 2) + std::pow(ys[i] - ys[i-1], 2)); + } + + CubicSpline phi(ts, xs); + CubicSpline psi(ts, ys); + + double xi; + for (int j = 0; j<=100; j++) { + xi = j*ts[12] / 100; + std::cout << phi.evaluate(xi) << " " << psi.evaluate(xi) << std::endl; + } +}