====== Lösungsvorschlag ======

==== Aufgabe 1 (Theorieaufgaben) ====
**a)**
  * O(n)
  * O(n)
  * O(n^2)
  * O(n^3)
  * O(n^2)
  * n-1
  * O(h^3)
  * O(h^4)

**b)**
z. B. Newton-Verfahren

**c)**
Anzahl korrekter Stellen nimmt pro Iterationsschritt quadratisch zu.

==== Aufgabe 2 (Direkte Verfahren fuer lin. Gleichungssysteme) ====
**a)**\\
                 | 0  1  0 |   |  4   6  4 |
  A = P_0 * A' = | 1  0  0 | * |  0   4  2 |
                 | 0  0  1 |   | -2  -1  3 |
      
           |   1  0 0 |   | 4  6  4 |
  = P_0 *  |   0  1 0 | * | 0  4  2 |
           | -1/2 0 1 |   | 0  2  5 |
      
           |   1    0   0 |   | 4  6  4 |
  = P_0 *  |   0    1   0 | * | 0  4  2 | = P_0 * L * R
           | -1/2  1/2  1 |   | 0  0  4 |
           
**b)**\\
                    |  0 |
                    |  1 |
  L * y = b <=> y = | -1 |
                    |  0 |
                             
                    |  1 |
                    |  0 |
  R * x = y <=> x = | -1 |
                    |  0 |
                    
 **c)**\\
            | -1  1 -1  1 |   | 2 0 -2 1 |
            |  1  1 -1 -1 |   | 0 2  1 0 |
  B = 1/2 * |  1  1  1  1 | * | 0 0  1 2 |
            |  1 -1 -1  1 |   | 0 0  0 1 |
            
==== Aufgabe 3 (Singulärwertzerlegung) ====
**a)**
  * V^T: Drehung von x
  * Sigma: Streckung / Stauchung von y
  * U: Drehung von z

**b)**
Die Singulärwerte sind die Quadratwurzeln der Eigenwerte von A^T*A bzw. A*A^T.

Sei A^T A = E_1 D_1 E_1^T und A A^T = E_2 D_2 E_2^T. Dann gilt: U = E_2, V = E_1.

**c)**
/A/_2 * /A^-1/_2 = 8 * 1 = 8 \\

**d)** 
  * im(A) = 1/5 * [(3, 0, 4, 0)^T, (0, -3, 0, -4)^T, (-4, 0, -3, 0)^T, (0, -4, 0, -3)^T]
  * ker(A) = {0} (der triviale Nullraum)
**e)**
            | 3 |                                  | 12 -6 -3 -6 |
            | 0 |                                  |  0  0  0  0 |
  8 * 1/5 * | 4 | * 1/5 * (4, -2, -1, -2) = 8/25 * | 16 -8 -4 -8 |
            | 0 |                                  |  0  0  0  0 |
                                             
==== Aufgabe 4 (Programmierung: Median Cut) ====
**a)**
  MedianCut::MedianCut(std::vector<Point2d>& pc):m_pointcloud(pc)
Oder:
  MedianCut::MedianCut(std::vector<Point2d>& pc){
    m_pointcloud = pc;
  }

**b)**
  void MedianCut::ComputeBoundingBox(const std::vector<Point2d>& pc, Point2d& min, Point2d& max) {
    min = new Point2d(pc[0]);
    max = new Point2d(pc[0]);
    
    for(auto &a : pc){
      if(a.x < min.x)
       min.x = a.x;
      if(a.y < min.y)
        min.y = a.y;
      if(a.x > max.x)
        max.x = a.x;
      if(a.y > max.y)
        max.y = a.y;
    }
  }

**c)**
  Point2d MedianCut::GetRepresentative(const std::vector<Point2d>& pointcloud) {
    Points2d result;
    for(auto &a : pointcloud){
      result += a;
    }
    
    result /= pointcloud.size();
    return result;
  }

**d)**
  std::vector<Point2d> MedianCut::ComputeMedianCut(int nCuts) {
    std::vector<Point2d> result;
    MedianCutRecursion(m_cloudpoint, 0, nCuts, result);
    return result;
  }

**e)**
  void MedianCut::MedianCutRecursion(std::vector<Point2d> subpointcloud, int cCuts, int nCuts, std::vector<Point2d>& result) {
    if (subpointcloud.empty()) return;
    std::vector<Point2d> sc = subpointcloud;
    
    if(sc.size() == 1 || cCuts >= nCuts) {
      result.push_back(GetRepresentative(sc));
      return;
    }
    
    Point2d min, max;
    ComputeBoundingBox(sc, min, max);
    
    if(max.x - min.x >= max.y - min.y)
      SortPointCloudAlongXAxis(sc);
    else
      SortPointCloudAlongYAxis(sc);
  
    std::size_t const half = sc.size() / 2;
    std::vector<Color> left(sc.begin(), sc.begin() + half);
    std::vector<Color> right(sc.begin() + half, sc.end());
  
    MedianCutRecursion(left, cCuts+1, nCuts, result);
    MedianCutRecursion(right, cCuts+1, nCuts, result);
  }

==== Aufgabe 5 (Iterative Loesungsverfahren) ====
**a)**\\
x_0 = (4 - (1 + 1)) / 2 = 1\\
x_1 = (2 - (1 + 1)) / 2 = 0\\
x_2 = (2 - (-1 - 1)) / 1 = 4\\
x_3 = (4 - (2 - 2)) / 2 = 2\\

**b)**\\
x_0 = (1 - 3/2) * 1 + (4 - (1 + 1)) / 2 = 1\\
x_1 = 1/2 + 3/2 * (2 - (1 + 1)) / 2 = 0.5\\
x_2 = -1/2 + 3/2 * (2 - (1/2 - 1)) / 1 = 3.25\\
x_3 = 1/2 + 3/2 * (4 - (2 + 1/2)) / 2 = 1.25\\

**c)**\\
  * Doppelte Gitterweite => 4-fache Anzahl Iterationen
  * Doppelte Fehlertoleranz => linear mehr Iterationen

   154 |  207 |  260
   620 |  830 | 1040
  2480 | 3320 | 4160
  
**d)**\\
  * Gauss-Seidel benoetigt halb so viele Iterationen wie Jacobi

    77 |  103 |  130
   310 |  415 |  520
  1240 | 1660 | 2080
  
 **e)**\\
  * SOR benoetigt Wurzel-n-mal so viele Iterationen wie Jacobi

  13 | 18 | 23
  25 | 35 | 45
  50 | 70 | 90

==== Aufgabe 6 (Konjugierte Richtungen, nichtlineare Optimierung) ====
**a)**\\
<note>Definition: zwei Vektoren u und v heißen A-konjugiert, wenn u^T * A * v = 0.</note>
Sei lambda eine reelle Zahl und sei r_1 = lambda *  (1, 0, 0)^T\\
r_1^T * A = (2, 1, 0) und (2, 1, 0) * (x, y, z)^T = 0 => r_2 = lambda * (1, -2, 0)^T\\
r_1^T * A = (2, 1, 0) und r_2^T * A = (0, -3, -2),\\
also (2, 1, 0) * (x, y, z)^T = 0 und (0, -3, -2) * (x, y, z)^T = 0 => r_3 = lambda *  (1/3, -2/3, 1)^T\\

**b)** Gradienten-Verfahren:\\
         | 2x - 2y - 4  |
  grad = | 2y^3 - 2x + 3 |

grad(3, 1) = (0, -1)^T

x_1 = x_0 + t*(-grad(3, 1)) = (3, 5/4)^T (wobei x_0 hier (x0,y0))

**c)** Newton-Verfahren:\\
      | 2      -2 |
  H = | -2   6y^2 |

            | 2   -2 |
  H(3, 1) = | -2   6 |
  
x_1 = x_0 - H(3, 1)^-1 * grad(3, 1) = 1/4 * (13, 5)^T

**d)**\\
  * max. Anzahl Iterationsschritte erreicht
  * gewuenschte Genauigkeit erreicht

==== Aufgabe 7 (Least Squares Approximation) ====
**a)**\\
[[https://de.wikipedia.org/wiki/Methode_der_kleinsten_Quadrate#Beispiel_mit_einer_Ausgleichsgeraden|Ansatz siehe Wikipedia]]\\
y = 3/5 * x + 2/5

**b)**\\
       |  4 |
  x = | -2 |
       | -2 |
       
==== Aufgabe 8 (Programmierung: Polynominterpolation) ====
**a)** \\
  // Wir bauen eine Vandermondematrix.
  
  uint32 c = controlPoints.size();
  Matrix m = new Matrix(c,c);
  for(int i = 0; i<c; i++){
    uint32 x = controlPoints[i].x;
    for(int k = 0; k<c; i++){
      m(i,k) = mul;
      mul*=mul;
    }
  }
  return m;
**b)** \\
  // pensi, preisi ist doch alles das selbe...

  uint32 c = controlPoints.size();
  Matrix b = new Matrix(c,1);
  for(int i = 0; i<c; i++){
    b(i,0) = controlPoints[i].y;
  }
  return b;
**c)** \\
  uint32 c = controlPoints.size();
  std::vector<float> n;
  
  // 4 -> Stammfunktion 4x^1, 4x -> 2x^2 also alter Wert/(i+1) wenn i bei 0 anfaengt
  
  for(int i = 0; i<c; i++){
    n.pushBack(f/i+1);
  }
  
  //Monom Kurve die die Ableitung repraesentiert
  
  MonomCurve tmp = new MonomCurve(n);
  return tmp.evaluateAt(b) - tmp.evaluateAt(a);
  
==== Aufgabe 9 (Bezier-Kurven) ====
**a)**\\
C_W(1/4) = (0, 1, 31/8)^T\\ 
C_O(1/4) = (4, 1, 23/8)^T\\

**b)**\\
S(s, t) = 1/2 * C_W =+ 1/2 * C_O = \\
= 1/2 * (b_0 * B(0, 2) + b_1 * B(1, 2) + b_2 * B(2, 2)) + 1/2 * (d_0 * B(0, 2) + d_1 * B(1, 2) + d_2 * B(2, 2)) =\\
= 1/2 * [ (b_0 + d_0) * B(0, 2) + (b_1 + d_1) * B(1, 2) + (b_2 + d_2) * B(2, 2) ] =\\
= (2,0,3)^T * B(0, 2) + (2,2,4)^T * B(1, 2) + (2,4,3)^T * B(2, 2)

**c)**\\
M(1/4) = (2, 1, 27/8)^T (z. B. <del>mithilfe De-Casteljau</del> viel einfacher: in a) hat man bereits C_W(1/4) und C_O(1/4) berechnet. Zu jedem t interpoliert S linear zwischen den Kurven! D. h. für s=1/2: S(1/2, 1/4) = 1/2 C_W(1/4) + 1/2 C_O(1/4))

**d)**\\
d_0 = S(0, 0) = C_W(0) = (0, 0, 4)\\
d_3 = S(1, 1) = C_O(1) = (4, 4, 4)

==== Aufgabe 10 (Multivariate Interpolation) ====

**10.1 a)**\\
Reihenfolge: links oben, rechts oben, links unten, rechts unten.
  * (rho > 0 and sigma > 0 and tau > 0) and (rho + sigma + tau = 1)
  * [(rho = 0 and 0 ≤  sigma, tau ≤  1) or (sigma = 0 and 0 ≤  rho, tau ≤  1) or (tau = 0 and 0 ≤  rho, sigma ≤  1)] and (rho + sigma + tau = 1)
  * rho = sigma
  * not (rho > 0 and sigma > 0 and tau > 0)

Beim zweiten hätte man auch jeweils eine der ≤ 1 Bedingungen weglassen können, denn diese folgen sowieso durch die Bedingung (rho + sigma + tau = 1).

**10.1 b)**\\
  * Q_0: (rho, sigma, tau) = (1/12, 1/4, 2/3)
  * Q_1: (rho, sigma, tau) = (1, 1, -1)

**10.1 c)**\\
  * a : b = beta : alpha
  * c : d = (1 - alpha) : alpha
  * e : f = gamma : alpha

**10.2**\\
  * P: in x-Richtung:
        * |P[00]P[10]| = |P[01]P[11]| = 6, |P[00]P| = |P[01]P| = 2, |P[11]P| = |P[10]P| = 4
        * f[x0] = 2/3 * 0 + 1/3 * 6 = 2, f[x1] = 2/3 * 6 + 1/3 * 2 = 14/3
        * 
        * in y-Richtung:
        * |P[x0]P[x1]| = 4, |P[x0]P| = 3, |P[x1]P| = 1
        * f = 3/4 * f[x1] + 1/4 * f[x0] = 4
  * M: 1/4 * 0 + 1/4 * 6 + 1/4 * 2 + 1/4 * 6 = 3,5