diff --git a/lecture.cls b/lecture.cls index 33a7c49..45c9d68 100644 --- a/lecture.cls +++ b/lecture.cls @@ -27,6 +27,17 @@ \RequirePackage{environ} \usetikzlibrary{quotes, angles} +\pgfplotsset{ + compat=1.15, + default 2d plot/.style={% + grid=both, + minor tick num=4, + grid style={line width=.1pt, draw=gray!10}, + major grid style={line width=.2pt,draw=gray!50}, + axis lines=middle, + enlargelimits={abs=0.2} + }, +} \DeclareOption*{\PassOptionsToClass{\CurrentOption}{article}} \DeclareOption{uebung}{ diff --git a/sose2020/num/uebungen/num9.pdf b/sose2020/num/uebungen/num9.pdf new file mode 100644 index 0000000..3da2932 Binary files /dev/null and b/sose2020/num/uebungen/num9.pdf differ diff --git a/sose2020/num/uebungen/num9.tex b/sose2020/num/uebungen/num9.tex new file mode 100644 index 0000000..6ad3c38 --- /dev/null +++ b/sose2020/num/uebungen/num9.tex @@ -0,0 +1,211 @@ +\documentclass[uebung]{../../../lecture} + +\title{Einführung in die Numerik: Übungsblatt 9} +\author{Leon Burgard, Christian Merten} + +\usepackage[]{subcaption} + +\begin{document} + +\punkte + +\begin{aufgabe} + \begin{enumerate}[a)] + \item Es ergibt sich durch Ausrechnen: + \begin{salign*} + a_0 &= y_0 = \frac{1}{2} \\ + a_1 &= \frac{y_1 - a_0 N_0(t_1)}{N_1(t_1)} = \frac{1 - \frac{1}{2}}{1 - \frac{1}{4}} = \frac{2}{3} \\ + a_2 &= \frac{y_2 - a_0 N_0(t_2) - a_1N_1(t_2)}{N_2(t_2)} = -\frac{4}{45} + \intertext{Damit folgt} + p_2(t) &= \frac{1}{2} + \frac{2}{3}\left( t - \frac{1}{4} \right) - \frac{4}{45} \left( t - \frac{1}{4} \right) (t- 1) + .\end{salign*} + \item Es folgt für $t_3 = 9$ $y_3 = 3$, also + \begin{salign*} + a_3 &= \frac{y_3 - a_0N_0(t_3) - a_1N_1(t_3) - a_2 N_2(t_3)}{N_3(t_3)} + = \frac{13}{1575} + \intertext{Damit folgt} + p_3(t) &= \frac{1}{2} + \frac{2}{3}\left( t - \frac{1}{4} \right) - \frac{4}{45} \left( t - \frac{1}{4} \right) (t- 1) + \frac{13}{1575} \left( t - \frac{1}{4} \right) (t-1)(t-4) + .\end{salign*} + \item Graph:\\ + \begin{tikzpicture} + \begin{axis}[default 2d plot, + xmin=0,xmax=10, + legend pos=outer north east] + \addplot[domain=0:10, smooth,green]{0.5 + 2/3*(x-1/4) - 4/45*(x-1/4)*(x-1)}; + \addlegendentry{$p_2(t)$} + \addplot[domain=0:10, smooth, blue]{0.5 + 2/3*(x-1/4) - 4/45*(x-1/4)*(x-1)+13/1575*(x-1/4)*(x-1)*(x-4)}; + \addlegendentry{$p_3(t)$} + \addplot[domain=0:10, smooth, red]{sqrt(x)}; + \addlegendentry{$\sqrt{t}$} + \end{axis} + \end{tikzpicture} + \end{enumerate} +\end{aufgabe} + +\begin{aufgabe} + \begin{enumerate}[a)] + \item Bezeichne: + \begin{align*} + &r_{i,0}(x) \coloneqq y_i \\ + &r_{i,k}(x) \coloneqq \frac{(x-x_i) p_{i+1, k-1}(x) - (x- x_{i+k})p_{i,k-1}(x)}{x_{i+k}-x_i} + .\end{align*} + Z.z.: $p_{i,0}(x) = r_{i,0}(x)$ und $p_{i,k}(x) = r_{i,k}(x)$. + + $r_{i,0}$ bzw. $r_{i,k}$ sind Polynome vom Grad $0$ bzw. $k$. D.h. es genügt zu zeigen, + dass sie die Interpolationseigenschaft erfüllen. Die Behauptung folgt dann aus der + Eindeutigkeit des Interpolationspolynoms. + + Für $r_{i,0}(x)$ ist die Interpolationseigenschaft trivialerweise für die eine Stützstelle + $x_i$ erfüllt, denn $r_{i,0}(x_i) = y_i$. + + Für $r_{i,k}$ gilt für $i \le j \le i+k$: + \begin{salign*} + r_{i,k}(x_{j}) &= \frac{(x_j-x_i) p_{i+1, k-1}(x_j) - (x_j- x_{i+k})p_{i,k-1}(x_j)}{x_{i+k}-x_i} \\ + &\stackrel{(*)}{=} \frac{(x_j - x_i) y_j - (x_j - x_{i+k})y_j}{x_{i+k}-x_i} \\ + &= y_j \frac{x_j - x_i - x_j + x_{i+k}}{x_{i+k}-x_i} \\ + &= y_j \frac{x_{i+k} - x_i}{x_{i+k}-x_i} \\ + &= y_j + .\end{salign*} + $(*)$: Falls $j = i$, dann ist $p_{i+1,k-1}(x_j)$ i.A. nicht $y_j$, aber dann ist + $(x_j - x_i) = (x_i - x_i) = 0$, also gilt dennoch $(x_j - x_i) p_{i+1,k-1}(x_j) = (x_j - x_i)y_j$. + Analog für $j = i+k$. + \item Durch Berechnung von $p_{0,3}(61.7)$ mit dem Schema aus (a) erhält man die Tageslänge + am Ort $E$ mit $19$h $29,55$m. + \end{enumerate} +\end{aufgabe} + +\begin{aufgabe} + \begin{enumerate}[a)] + \item + \begin{itemize} + \item Für die Koeffizienten $a_i$ gilt $a_i = y_i$. Also keine Operationen nötig. + \item Auswertung von $L_i^{(n)}(\xi) = \prod_{j=0,j\neq i}^{n} \frac{(\xi-x_j)}{x_i - x_j} $: + $n+1$ Faktoren mit $3$ Operationen plus $n$ Multiplikationen für das Produkt. Gesamt: + $3(n+1) + n = 4n+4+n=5n+4=5(n+1)-1$. + + Auswertung von $p(\xi)$: $n+1$ Summanden mit $5(n+1)-1$ Operationen plus Multiplikation + mit $y_i$, ergibt $(n+1)\cdot (5(n+1)-1+1) = 5(n+1)^2$. Mit zusätzlich $n$ Additionen + für die Auswertung der Summe ergibt sich insgesamt $5(n+1)^2 + n = \mathcal{O}(n^2)$. + \end{itemize} + \item + \begin{itemize} + \item + Berechnung der Koeffizienten erfordert die Lösung eines LGS mit unterer Dreicksmatrix. + Dies erfordert $\mathcal{O}^2$ Operationen. + \item Auswertung von $N_i(\xi) = \prod_{j=0}^{i-1} (\xi - x_j) $ erfordert $1$ Addition pro Faktor und + insgesamt $i-1$ Multiplikationen für das Produkt, also insgesamt: $2i - 1$. + + Auswertung von $p(\xi)$ erfordert entsprechend die Auswertung von $N_i(\xi)$ und die + Multiplikation mit dem Koeffizienten und Summation über alle $(n+1)$ Summanden. Ergibt also + insgesamt mit kleinem Gauß $(n+1)(n+2) = n^2 + 3n + 2 = \mathcal{O}(n^2)$. + \end{itemize} + \item + \begin{itemize} + \item Berechnung der Koeffizienten erfordert die Lösung eines vollen LGS, also + $\mathcal{O}(n^{3})$ Operationen. + \item Auswertung eines Summanden: $a_i t ^{i}$ erfordert $i-1 + 1 = i$ Multiplikationen. + + Auswertung von $p(x)$ erfordert das Summieren von $n+1$ Summanden, also $n$ zusätzliche Additionen. + Insgesamt folgt also wieder mit kleinem Gauß $\frac{(n+1)(n+2)}{2} + n = \frac{n^2 + 3n + 2}{2} + n = \frac{n^2 + 5n +4}{2} = \mathcal{O}(n^2)$. + \end{itemize} + \item Auswertung von $p_{i,k}$ im Neville Schema erfordert + \begin{align*} + N(k) &= 2 + N(k-1) + 1 + 1 + 1 + N(k-1) + 1 + 1 \\ + &= 7 + 2 N(k-1) \\ + &= 7 + 2 \cdot 7 + 4 \cdot N(k-2) = \ldots = \sum_{j=0}^{k-1} 7 \cdot 2^{j} \\ + &= 7 \frac{2^{k} - 1}{2 - 1} \\ + &= 7 (2^{k} - 1) \\ + &= \mathcal{O}(2^{k}) + .\end{align*} + \end{enumerate} +\end{aufgabe} + +\begin{aufgabe} + \begin{enumerate}[a)] + \item Auszug aus \textit{polynom.cc} + \begin{lstlisting}[language=C++, title=Auswertung des Interpolationspolynoms für vorgegebene Stüztstellen, captionpos=b] +// Auswertung von p_{i,k}(t) mit Neville Schema bei t=x +template +REAL evaluateNeville(REAL x, std::vector &xs, std::vector &ys, int i, int k) { + if(k == 0) { + return ys[i]; + } else { + // nevile rekursionsformel + return ((x - xs[i])*evaluateNeville(x, xs, ys,i+1,k-1) - (x - xs[i+k])*evaluateNeville(x, xs, ys, i, k-1))/(xs[i+k]-xs[i]); + } +} + +// Auswertung eines Interpolationspolynoms p(t) bei t=x +template +REAL evaluate(REAL x, std::vector &xs, std::vector &ys) { + // nutze neville verfahren mit p(x) = p_{0,n}(x) + evaluateNeville(x, xs, ys, 0, xs.size()-1); +}\end{lstlisting} +\item Auszug aus \textit{polynom.cc} +\begin{lstlisting}[language=C++, title=Auswertung von I-Polynomen für äquidistante Stützstellen, captionpos=b] +// Auswertung eines Interpolationspolynoms einer Funktion bei vorgegebenen Stuetzstellen +template +REAL evaluateFunction(REAL x, std::vector &xs, REAL(*f)(REAL)) { + std::vector ys(xs.size()); + for (int i = 0; i < xs.size(); i++) { + // berechne stuetzstellen + ys[i] = f(xs[i]); + } + // verwende polynominterpolation + return evaluate(x, xs, ys); +} + +// Auswertung eines Interpolationspolynoms einer Funktion bei aequidistanten Stuetzstellen +template +REAL evaluateFunctionEqualDist(REAL x, REAL a, REAL b, REAL h, REAL(*f)(REAL)) { + std::vector xs(std::floor((b-a)/h)); + std::vector ys(xs.size()); + for (int i = 0; i < xs.size(); i++) { + // berechne stuetzstellen + xs[i] = a + i*h; + // und funktionswerte + ys[i] = f(xs[i]); + } + // werte polynom aus + return evaluate(x, xs, ys); +}\end{lstlisting} +\begin{figure}[h] + \begin{subfigure}[b]{.5\linewidth} + \begin{tikzpicture} + \begin{axis}[default 2d plot, xmin=-1, xmax=1] + \addplot[green] table {f1_5.dat}; + \addlegendentry{$n=5$} + \addplot[orange] table {f1_10.dat}; + \addlegendentry{$n=10$} + \addplot[blue] table {f1_20.dat}; + \addlegendentry{$n=20$} + \addplot[red, samples=100] {1/(1+x^2)}; + \addlegendentry{$\frac{1}{1+x^2}$} + \end{axis} + \end{tikzpicture} + \caption{$f_1 = \frac{1}{1+x^2}$} + \end{subfigure} + \begin{subfigure}[b]{.5\linewidth} + \begin{tikzpicture} + \begin{axis}[default 2d plot, xmin=-0.5, xmax=0.5, legend pos=outer north east, + restrict y to domain=-1:5] + \addplot[green] table {f2_5.dat}; + \addlegendentry{$n=5$} + \addplot[orange] table {f2_10.dat}; + \addlegendentry{$n=10$} + \addplot[blue] table {f2_20.dat}; + \addlegendentry{$n=20$} + \addplot[red, samples=5000] {sqrt(abs(x))}; + \addlegendentry{$\sqrt{|x|} $} + \end{axis} + \end{tikzpicture} + \caption{$f_2 = \sqrt{|x|} $} + \end{subfigure} +\end{figure} +\item Für $\frac{1}{1+x^2}$ funktioniert die Interpolation mit äquidistanten Stützstellen sehr gut. + Für $\sqrt{|x|} $ entstehen durch die nicht differenzierbare Stelle bei $x=0$ sehr große Abweichungen, + die zu den Rändern mit wachsendem Polynomgrad sogar zunehmen. + \end{enumerate} +\end{aufgabe} + +\end{document} diff --git a/sose2020/num/uebungen/polynom.cc b/sose2020/num/uebungen/polynom.cc new file mode 100644 index 0000000..6bbc407 --- /dev/null +++ b/sose2020/num/uebungen/polynom.cc @@ -0,0 +1,75 @@ +#include +#include +#include + +// Auswertung von p_{i,k}(t) mit Neville Schema bei t=x +template +REAL evaluateNeville(REAL x, std::vector &xs, std::vector &ys, int i, int k) { + if(k == 0) { + return ys[i]; + } else { + // nevile rekursionsformel + return ((x - xs[i])*evaluateNeville(x, xs, ys,i+1,k-1) - (x - xs[i+k])*evaluateNeville(x, xs, ys, i, k-1))/(xs[i+k]-xs[i]); + } +} + +// Auswertung eines Interpolationspolynoms p(t) bei t=x +template +REAL evaluate(REAL x, std::vector &xs, std::vector &ys) { + // nutze neville verfahren mit p(x) = p_{0,n}(x) + evaluateNeville(x, xs, ys, 0, xs.size()-1); +} + +// Auswertung eines Interpolationspolynoms einer Funktion bei vorgegebenen Stuetzstellen +template +REAL evaluateFunction(REAL x, std::vector &xs, REAL(*f)(REAL)) { + std::vector ys(xs.size()); + for (int i = 0; i < xs.size(); i++) { + // berechne stuetzstellen + ys[i] = f(xs[i]); + } + // verwende polynominterpolation + return evaluate(x, xs, ys); +} + +// Auswertung eines Interpolationspolynoms einer Funktion bei äquidistanten Stützstellen +template +REAL evaluateFunctionEqualDist(REAL x, REAL a, REAL b, REAL h, REAL(*f)(REAL)) { + std::vector xs(std::floor((b-a)/h)); + std::vector ys(xs.size()); + for (int i = 0; i < xs.size(); i++) { + // berechne stuetzstellen + xs[i] = a + i*h; + // und funktionswerte + ys[i] = f(xs[i]); + } + // werte polynom aus + return evaluate(x, xs, ys); +} + +// beispiel funktionen +template +REAL f1(REAL x) { + return 1/(1+std::pow(x,2)); +} + +template +REAL f2(REAL x) { + return std::sqrt(std::abs(x)); +} + +int main() { + // polynom grad + int n = 5; + // abstand zwischen stützstellen + double h = 2.0 / n; + + // schrittweite + double dx = 2.0 / 1000; + double x = 0; + for(int i = 0; i < 1000; i++) { + x = -1.0 + i*dx; + // werte polynom mit äquidistanten stuetzstellen zwischen -1 und 1 mit abstand h bei t=x aus + std::cout << x << " " << evaluateFunctionEqualDist(x, -1.0, 1.0, h, f2) << std::endl; + } +}