diff --git a/sose2020/num/uebungen/iterative_plot.plt b/sose2020/num/uebungen/iterative_plot.plt new file mode 100644 index 0000000..f3e04d1 --- /dev/null +++ b/sose2020/num/uebungen/iterative_plot.plt @@ -0,0 +1,19 @@ +set terminal png size 1000,1000 +set output 'iterative_solvers_plot.png' + +set logscale xy + +set ylabel "Laufzeit in Sekunden" +set xlabel "Matrixgröße" + +set format y "%g" + +plot 'richardson_1.dat' title "Richardson Iteration v1", \ + 'richardson_2.dat' title "Richardson Iteration v2", \ + 'jacobi_1.dat' title "Jacobi Iteration v1", \ + 'jacobi_2.dat' title "Jacobi Iteration v2",\ + 'gauss_seidel_1.dat' title "Gauss Seidel Iteration v1",\ + 'gauss_seidel_2.dat' title "Gauss Seidel Iteration v2",\ + 'linsolve_1.dat' title "hdnum::linsolve",\ + 0.000000020*x**4 title "c n^4", \ + 0.000000025*x**3 title "b n^3" diff --git a/sose2020/num/uebungen/iterative_solvers_plot.png b/sose2020/num/uebungen/iterative_solvers_plot.png new file mode 100644 index 0000000..78a8945 Binary files /dev/null and b/sose2020/num/uebungen/iterative_solvers_plot.png differ diff --git a/sose2020/num/uebungen/num7.pdf b/sose2020/num/uebungen/num7.pdf new file mode 100644 index 0000000..28552d2 Binary files /dev/null and b/sose2020/num/uebungen/num7.pdf differ diff --git a/sose2020/num/uebungen/num7.tex b/sose2020/num/uebungen/num7.tex new file mode 100644 index 0000000..009e60c --- /dev/null +++ b/sose2020/num/uebungen/num7.tex @@ -0,0 +1,175 @@ +\documentclass[uebung]{../../../lecture} + +\title{Einführung in die Numerik: Übungsblatt 7} +\author{Leon Burgard, Christian Merten} + +\begin{document} + +\punkte + +\begin{aufgabe}[] + \begin{enumerate}[a)] + \item Mit der Block LU-Zerlegung von $A$ folgt + \begin{align*} + A &= \begin{pmatrix} Id & 0 \\ + A_{21}A_{11}^{-1} & Id \end{pmatrix} + \begin{pmatrix} A_{11} & A_{12} \\ + 0 & S\end{pmatrix} + = \begin{pmatrix} + A_{11} & A_{12} \\ + A_{21} A_{11}A_{11}^{-1} A_{22}A_{11}^{-1}A_{12} + S + \end{pmatrix} + = \begin{pmatrix} + A_{11} & A_{12} \\ + A_{21} & A_{21}A_{11}^{-1}A_{12} + S + \end{pmatrix} + .\end{align*} + Damit folgt + \[ + A_{22} = A_{21}A_{11}^{-1}A_{12} + S \implies S = A_{22} - A_{21}A_{11}^{-1}A_{12} + .\] + \item Sei $A$ hermitesch und positiv definit. Dann ist + \[ + \overline{A}^{T} = \begin{pmatrix} \overline{A_{11}}^{T} & \overline{A_{21}}^{T} \\ + \overline{A_{12}}^{T} & \overline{A_{22}}^{T} + \end{pmatrix} = + \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} = A + .\] + Damit folgt $\overline{A_{11}}^{T} = A_{11}$ und $\overline{A_{12}}^{T} = A_{21}$. Dann folgt + \[ + \overline{S}^{T} = \overline{A_{22}}^{T} - \overline{A_{12}}^{T} \overline{A_{11}}^{-T} \overline{A_{21}}^{T} + = A_{22} - A_{21}A_{11}^{-1}A_{12} = S + .\] Also $S$ und $A_{11}$ hermitesch. Da $A$ hermitesch und positiv definit, sind alle + führenden Hauptminoren positiv, d.h. auch alle führenden Hauptminoren von $A_{11}$ sind positiv, + d.h. $A_{11}$ positiv definit. + \end{enumerate} +\end{aufgabe} + +\begin{aufgabe} + Beh.: Der Algorithmus ist wie angegeben durchführbar. + \begin{proof} + Induktionsbehauptung: $\forall 1 \le j < n$: $u_j \neq 0$ und $|u_j| > |b_j|$. + Endliche Induktion über $j < n$. $j = 1$: $u_1 = a_1 \neq 0$. $|u_1| > |b_1|$. + + Sei $j < n$ und Induktionsbehauptung für $j-1$ gezeigt. Dann gilt $u_{j-1} \neq 0$, also + $l_j = \frac{c_j}{u_{j-1}}$ und $u_j = a_j - \frac{c_j}{u_{j-1}}b_{j-1}$. Damit folgt + \begin{salign*} + |u_j| &= \left| a_j - \frac{c_j}{u_{j-1}}b_{j-1} \right| \\ + &\ge \left| |a_j| - \frac{|c_j|}{|u_{j-1}|}|b_{j-1}| \right| \\ + &\ge \Big| |b_j| + \underbrace{|c_j|}_{\neq 0} \Big( 1 - \underbrace{\frac{|b_{j-1|}}{|u_{j-1}|}}_{ \text{I.V.:} < 1} \Big) \Big| \\ + &> |b_j| > 0 + .\end{salign*} + Das zeigt die Induktionsbehauptung. + + Für $j = n$ folgt ganz analog + \begin{salign*} + |u_n| &= \left| a_n - \frac{c_n}{u_{n-1}}b_{n-1} \right| \\ + &\ge \Big| |c_n| \Big(1 - \underbrace{\frac{|b_{n-1}|}{|u_{n-1}|}}_{< 1}\Big) \Big| \\ + &> 0 + .\end{salign*} + \end{proof} + Beh.: Der Algorithmus liefert die angegebene LU-Zerlegung. + \begin{proof} + Es müssen je Zeile nur die $c_j$ eliminiert werden. Das wird mit $l_j = c_j / u_{j-1}$ erreicht. + Da $b_j \neq 0$, wird noch die Diagonale modifiziert um $-l_{j}b_{j-1}$. Die $b_j$ werden + nicht verändert, da die Elemente in $A$ oberhalb der $b_j$ null sind. Damit folgt die angegebene + LU-Zerlegung. + \end{proof} + Beh.: $\text{det}(A) \neq 0$. + \begin{proof} + Es ist $\text{det}(A) = \text{det}(L) \cdot \text{det}(U) = 1 \underbrace{\cdot u_1 \cdots u_n}_{\neq 0} \neq 0$ + \end{proof} +\end{aufgabe} + +\begin{aufgabe} + Sei $A \in \R^{n \times n}$ gegeben durch + \[ + a_{ij} = \begin{cases} + +1 & i=j \lor j=n \\ + -1 & i > j \\ + 0 &\text{sonst} + \end{cases} + .\] + \begin{enumerate}[a)] + \item Es gilt nach VL für $1 \le k < n$: + \[ + l_i^{(k)} = \begin{cases} + 0 & 1 \le i \le k \\ + a_{ik}/a_{kk} & k < i \le n + \end{cases} + .\] + Für $1 \le i < n$ gilt: Wegen $a_{ij} = 0$ für $j > i$, werden die Pivotelemenete $a_{ii} = 1$ + nicht modifiziert. Damit folgt + \[ + l_{ij} = \begin{cases} + 0 & i < j \\ + +1 & i = j \\ + -1 & i> j + \end{cases} + ,\] also $|l_{ij}| \le 1$. Da $l_{ij} = -1$ für $i > j$, gilt für $i > k$: + \[ + a_{in}^{(k)} = a_{in}^{(k-1)} - (-1) \cdot a_{kn}^{(k)} + .\] Damit folgt + \[ + u_{in}^{(k)} = \begin{cases} + u_{in}^{(k-1)} & 1 \le i \le k \\ + u_{in}^{(k-1)} + u_{kn}^{(k)} & k+1 \le i \le n + \end{cases} + = \begin{cases} + u_{in}^{(k-1)} & 1 \le i \le k \\ + 2 u_{in}^{(k-1)} & k+1 \le i \le n + \end{cases} + .\] Insgesamt folgt + \[ + u_{nn} = u_{n n}^{(n-1)} = 2 u_{nn}^{(n-2)} = \ldots = 2^{n-1} u_{nn}^{(1)} = 2^{n-1} + .\] + \item Verwende Spaltenvertauschungen $Q \in \R^{n \times n}$, $j$-te Spalte von $Q$ gegeben als + \[ + Q_j = \begin{cases} + e_n & j = 1 \\ + e_{j-1} & 1 < j \le n + \end{cases} + .\] Dann hat $AQ$ die Form + \[ + (AQ)_{ij} = \begin{cases} + +1 & j = 1 \lor i+1 = j \\ + -1 & j\neq 1 \land j < i+1 \\ + 0 & \text{sonst} + \end{cases} + .\] + Induktion über $n$. Für $n = 2$ gilt + \[ + A = \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} + \implies + AQ = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} + = + \underbrace{\begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}}_{=: L} + \underbrace{\begin{pmatrix} 1 & 1 \\ 0 & -2 \end{pmatrix} }_{=: U} + .\] Es gilt also $u_{nn} = -2$, insbes. $|u_{nn}| = 2 = \max \{|1|, |-2|\} $. + + Sei nun $n \in \N$ beliebig und Beh. gezeigt für $n-1$. Dann ist + \[ + AQ = \begin{bmatrix} \tilde{A} & \begin{matrix} 0 \\ \vdots \\ 1 \end{matrix} \\ + \begin{matrix} 1 & -1 & \cdots & -1 \end{matrix} & -1 \end{bmatrix} + .\] Wende Operationen der Gauß-Elimination von $\tilde{A} \in \R^{(n-1)\times (n-1)}$ auf $AQ$ an + Für $1 \le i \le n-2$ gilt $a_{in}= 0$, also bleibt $n$-te Spalte unverändert. Letzte Zeile + zu $0$ eliminiert, bis auf $a_{n(n-1} = -1$. Nach I.V. gilt + jetzt $a'_{(n-1)(n-1)} = -2$. Mit $l = 1$ folgt + \[ + u_{nn} = a'_{nn} = a_{nn} - a_{(n-1)n} = -1 -1 = -2 + .\] Es ist weiter + \[ + u_{ni} = a_{ni} = \begin{cases} + 0 & 1 \le i < n -1 \\ + 1 & i = n-1 \\ + -2 & i = n + \end{cases} + .\] Damit folgt die Behauptung. + \end{enumerate} +\end{aufgabe} + +\begin{aufgabe} + siehe \textit{prog\_iterative\_solvers.cc} und \textit{iterative\_solvers\_plot.png}. +\end{aufgabe} + +\end{document} diff --git a/sose2020/num/uebungen/prog_iterative_solvers.cc b/sose2020/num/uebungen/prog_iterative_solvers.cc new file mode 100644 index 0000000..99c87f1 --- /dev/null +++ b/sose2020/num/uebungen/prog_iterative_solvers.cc @@ -0,0 +1,183 @@ +#include // notwendig zur Ausgabe +#include +#include "hdnum.hh" // hdnum header + +namespace hdnum { + + template + class SparseMatrix { + + struct MatrixEntry { + int i; + int j; + REAL value; + }; + + public: + + void AddEntry (int i, int j, REAL value) { + assert(i >= 0); + assert(j >= 0); + if (value != .0) + entries.push_back(MatrixEntry{.i=i, .j=j, .value=value}); + } + + template + void mv (Vector& y, const Vector& x) { + + zero(y); + + for (MatrixEntry& matrix_entry : entries) { + assert(y.size() > matrix_entry.i); + assert(x.size() > matrix_entry.j); + y[matrix_entry.i] += matrix_entry.value * x[matrix_entry.j]; + } + } + + //private: + std::vector entries; + }; + +} + +// generic iterative updater function +typedef void iterativeUpdater(hdnum::DenseMatrix& A, + hdnum::SparseMatrix& A_s, + int N, + hdnum::Vector& x_tmp, + hdnum::Vector& b, + hdnum::Vector& d, + hdnum::Vector& x); + +// solve iterative based on given method +void solveIterative(hdnum::DenseMatrix& A, + hdnum::SparseMatrix& A_s, + hdnum::Vector & x, + hdnum::Vector& b, + iterativeUpdater update) { + int N = A.rowsize(); + // defekt vektor + hdnum::Vector d(N); + // temporary variable + hdnum::Vector x_tmp(N); + A_s.mv(x_tmp, x); + d = b - x_tmp; + double initialDefekt = d.two_norm(); + while (d.two_norm() > initialDefekt * 10e-4) { + // run one iteration + update(A, A_s, N, x_tmp, b, d, x); + } +} + +// richardson: W = 1/omega I => W^-1 = omega I +inline void richardson(hdnum::DenseMatrix& A, + hdnum::SparseMatrix& A_s, + int N, + hdnum::Vector& x_tmp, + hdnum::Vector& b, + hdnum::Vector& d, + hdnum::Vector& x) { + A_s.mv(x_tmp, x); + d = b - x_tmp; + // x(k+1) = x(k) + omega * d(k) + // hier omega = -0.01 + for (int i=0; i& A, + hdnum::SparseMatrix& A_s, + int N, + hdnum::Vector& x_tmp, + hdnum::Vector& b, + hdnum::Vector& d, + hdnum::Vector& x) { + A_s.mv(x_tmp, x); + d = b - x_tmp; + // jacobi iteration, teile durch Diagonal Elemente + // x(k+1) = x(k) + D^-1 d(k) + for (int i=0; i& A, + hdnum::SparseMatrix& A_s, + int N, + hdnum::Vector& x_tmp, + hdnum::Vector& b, + hdnum::Vector& d, + hdnum::Vector & x) { + A_s.mv(x_tmp, x); + d = b - x_tmp; + x_tmp = d; + // x(k+1) = x(k) + v(k) + // W*v(k) = d(k) + // W untere Dreicksmatrix, loese W*v(k) = d(k) durch Vorwaertseinsetzen + for (int i = 0; i& A, + hdnum::SparseMatrix& A_s, + hdnum::Vector& x, + hdnum::Vector& b, + iterativeUpdater updater) { + int N = A.rowsize(); + hdnum::Vector y(N); + // teste iteration mit gegebenem verfahren + solveIterative(A,A_s,x,b,updater); + // vergleiche mit LU Zerleger + hdnum::linsolve(A,y,b); + std::cout << "Relativer Fehler: " << (y - x).two_norm() / y.two_norm() << std::endl; +} + +int main () { + for(int n=4; n<=9; n++) { + int N = pow(2,n); + + // Testmatrix aufsetzen + hdnum::DenseMatrix A(N,N,.0); + hdnum::SparseMatrix A_s; + for (typename hdnum::DenseMatrix::size_type i=0; i 0) { + A[i][i-1] = 1.0; + A_s.AddEntry(i, i-1, 1.0); + } + if (i + 1 < A.colsize()) { + A[i][i+1] = 1.0; + A_s.AddEntry(i, i+1, 1.0); + } + A[i][i] -= 2.0; + A_s.AddEntry(i, i, -2.0); + } + + // Rechte Seite und Lösungsvektor + hdnum::Vector x(N, 0.0); + hdnum::Vector b(N, 1.0); + + // Lösen Sie nun A*x=b iterativ + + // Pretty-printing fuer Vektoren + x.scientific(false); + x.width(15); + + // testApproximation(A,A_s,x,b,richardson); + + hdnum::Timer myTimer; + solveIterative(A,A_s,x,b,gaussSeidel); + //hdnum::linsolve(A,x,b); + std::cout << N << " " << myTimer.elapsed() << std::endl; + } +}