kit-books

summary.typ (23798B)

---

 #set page("a4", margin: 1cm, numbering: "1/1", header: rect(width: 100%, stroke:
 (bottom: 1pt), [Numerik Zusammenfassung #h(1fr) Gero Beckmann]))
 
 #let TODO(..rest, content) = text(stroke: red, ..rest, [TODO: #content])
 
 == Gleitpunktzahlen
 
 $
 "FL" = {
   +- sum_(l=1)^(L_m) a_l B^(-l) B^e : e = e_"min" + sum_(l=0)^(L_e - 1) c_l B^l, 
   a_l, c_l in {0, ..., B-1}, a_1 != 0
 } union {0} 
 $
 
 #columns()[
 $
 & e_"max" = e_"min" + B^(L_e) - 1 \
 & max "FL" = - min "FL" = B^(e_max) ( 1 - B^(-L_m)) \
 & min "FL"_+ = B^(e_min - 1) \
 & B^(-1) <= abs(m) < 1
 $
 
 #colbreak()
 
 Relative Maschinengenauigkeit 
 
 $
 abs(x - "fl"(x)) / abs(x) <= (B^(1- L_m)) / 2 =: "eps" \
 "fl"(x) = x ( 1 + epsilon) " mit " abs(epsilon) <= "eps" \
 x hat(compose) y := "fl"(x compose y) = (x compose y)(1 + epsilon) " mit " abs(epsilon) <= "eps"
 $
 ]
 
 Stabilität eines Algorithmus: Fehlerverstärkung des Verfahrens ist "moderat"
 groß (im vergleich zum unvermeidbaren Fehler). 
 
 #columns(gutter: 0.1cm)[
 == Cholesky-Zerlegung
 
 Ist eine $A in bb(R)^(N times N)$ symmetrisch und positix definit, existirt eine
 Cholesky-Zerlegung $A = L L^T$
 
 1. Berechne Cholesky-Zerlegung $A = L L^T$
 2. Löse $L y = b$ durch Vorwärtssubstitution
 3. Löste $L^T x = y$ durch Rückwärtssubstitution
 
 $l_(n n) > 0$ macht die Zerlegugn eindeutig. 
 
 $
 Q(lambda m) = lambda Q m = - lambda w " für alle " lambda in bb(R) \
 Q y = y " für alle " y in bb(R)^N " mit " w^T y = 0 \
 "Spiegelebene " E = {y in bb(R)^M: w^T y = 0}
 $
 #colbreak()
 == LR-Zerlegung
 
 Ist jeder Hauptminor von $A in bb(R)^(N times N)$ regulär, existire eine
 LR-Zerlegen $A = L R$.  Ist $A in bb(R)^(N times N)$ regulär existiert eine LR-Zerlegung mit
 Spaltenpivotwahl $P A = L R$. 
 
 1. Berechne Zerlegung $P A = L R$ durch Gauß-Elimination
 2. Löse $L y = P b$ durch Vorwärtssubstitution in $1/2 N^2$ Ops
 3. Löse $R x = y$ durch Rückwärtssubstitution $x_n = [y_n - sum_(j=n+1)^N r_(n j) x_j] / r_(n n)$ in $O(1/2 N^2)$ Ops
 
 $l_(n n) = 1$ liefiert Eindeutigkeid der Zerlegung. 
 ]
 
 #pagebreak()
 
 == Kondition und Normen
 
 Problem: Wie wirken sich Störungen der Eingabegrößen (hier A und b) auf die
 Lösung unabhängig vom gewählten Algorithmus aus?
 
 === Zugehörige Matrixnorm
 
 $
 norm(A) &:= sup_(x != 0) norm(A x)/norm(x) = sup_(x !=0) norm(A x/norm(x)) =
 sup_(norm(x)=1) norm(A x) = max_(norm(x) = 1) norm(A x) \
 
 norm(A^(-1)) &= sup_(x != 0) norm(A^(-1) x)/ norm(x) =^(A^(-1)x = 1) sup_(z !=
 0) norm(z)/norm(A z) = max_(norm(z) = 1) 1/norm(A z) = 1 / (min_(norm(z) = 1)
    norm(A z))
 $
 
    #table(columns: 4, align: center + horizon,
 [Spaltensummennorm], 
 $ norm(A)_1 = max_(m=1,...,N) sum_(n=1)^N abs(a_"nm") $,
 $ norm(x)_1 = sum_(n=1)^N abs(x_n) $, 
 [1-Norm],
 [Zeilensummennorm], 
 $ norm(A)_oo = max_(n=1,...,N) sum_(m=1)^N abs(a_"nm") $,
 $ norm(x)_oo = max_(n=1,...,N) abs(x_n) $, [Maximumsnorm],
 [Spektralnorm], $ norm(A)_2 = sqrt("größter EW von " A^T A) $, $ norm(x)_2 = sqrt(sum_(n=1)^N x_n^2) $, [euklidische Norm]
  )
 
  Eigenschaften: positive Definitheit $norm(x) >= 0$, Dreiecksungleichung $norm(x
  + y ) <= norm(x) + norm(y)$, Homogenität $norm(lambda x) = abs(lambda) norm(x)$
 
 $norm(A x) <= norm(A) norm(x)$, $norm(I_N) = 1$, Submiltiplikativität $norm(A B)
 <= norm(A) norm(B)$, $norm(x) = sqrt(angle.l x"," x angle.r)$
 
 === Störungen
 
 #let cond = "cond";
 
 $
 x - tilde(x) &= A^(-1) b - A^(-1) tilde(b) = A^(-1)(b - tilde(b)) \
 <=> norm(x - tilde(x)) &= norm(A^(-1) (b - tilde(b))) <= norm(A^(-1)) norm(b - tilde(b))
 && "(absolute Abweichung)" \
 <=> norm(x - tilde(x)) / norm(x) &= norm(A^(-1) ( b - tilde(b))) / norm(x) <=
 (norm(b) norm(A^(-1))) / norm(x) norm(b - tilde(b)) / norm(b) && "(relative Abweichung)" \
 <=> norm(x - tilde(x)) / norm(x) &<=^(norm(b) = norm(A x) <= norm(A) norm(x))
 underbrace(norm(A) norm(A^(-1)), =: cond(A)) norm(b - tilde(b)) / norm(b) $
 
 === Konditionszahl einer Matrix
 
 $
 #rect($cond(A) := norm(A) norm(A^(-1))$) >=^"Submiltiplikativit" norm(A A^(-1))
 = norm(I_N) = 1
 $
 
 $
 cond_2(A) = (max{abs(lambda) : lambda "EW von " A})/(max{abs(lambda) : lambda
 "EW von " A}), cond(A) = (max_(norm(y) = 1) norm(A y))/(min_(norm(z) = 1) norm(
   A z)), cond(A) = cond(alpha A), alpha in bb(R) "\\" {0}
 $
 
 Seien relative Abweichung beschränkt
 $
 norm(A - tilde(A)) / norm(A) <= epsilon_A, "       " && norm(b- tilde(b)) / norm(b) <= epsilon_b
 $
 Dann gilt im Fall $epsilon_A cond(A) < 1$ die Abschätzung
 $
 norm(x - tilde(x))/norm(x) <= cond(A) / (1 - epsilon_A cond(A)) (epsilon_A + epsilon_b)
 $
 
 $cond(alpha I_N) = 1$, für orthogonale Matrizen Q: $ norm(Q x)_2^2 = x^T Q^T T
 x = x^T x = norm(x)_2^2, norm(Q)_2 = 1 => cond_2(Q) = 10 = norm(Q) norm(Q^(-1)) =
 norm(Q)^2$
 
 #pagebreak()
 
 == QR-Zerlegung
 
 Geg. $A in bb(R)^(M times N), M >= N, "voller Rang" (=N)$\
 Ges. orthogonales $Q in bb(R)^(M times N)$ und $R = mat(dash(R);0) in bb(R)^(M
 times N)$, rechte obere Dreiecksmatrix $dash(R) in bb(R)^(N times N)$ 
 
 #rect()[
 === Householder-Transformation
 
 Householdermatrix Q ist symmetrisch, orthogonal und eine Spiegelung eines
 Vektors v auf ein Vielfaches des Vektors u $(=e^1)$ mit Spiegelachse $perp w$ (Householdervektor). 
 
 $
 Q = I_N 2 w w^T => Q v = sigma e^1 = -"sign"(v_1) norm(v)_2 e^1 = v - 2 w^T v w
 #h(1cm) (2M " Ops")
 $
 
 Q wird eigl. nie ausgerechnet, sondern nur w gespeichert und QA wird berechnet
 aus $(Q a^1, Q a^2, ...)$ Spaltenvektoren von A. (2M(N-1)) Ops
 
 $
 "(Householdervektor) " w = (v - u)/norm(v - u)_2 = (v - sigma e^1)/norm(v - sigma e^1)_2 
 #h(1cm) sigma = -"sign"(v_1)norm(v)_2 = cases(-&norm(v)_2 ", falls " v_1 > 0,
 &norm(v)_2 ", falls " v_1 <= 0)
 $
 ]
 
 Berechne Householder-Transformation des ersten Spaltenvektors $a^1$ von $A^((n))$ 
 
 $
 w = (a^1 - sigma e^1) / norm(a^1 - sigma e^1) 
 #h(1cm) , 
 Q a^1 = r^1 = sigma e^1 = -"sign"(a^1) norm(a^1)_2 e^1
 #h(1cm) , 
 Q^((n)) A^((n)) 
 = mat(
 r_(11), r_(12)^T;
 , A^((n+1))
 )
 = mat(
 sigma, r_(12)^T;
 0, A^((n+1))
 )
 $
 
 Durch die Householder-Transformation bekommen wir die Matrix für die nächste
 Iteration. 
 
 Householdervektor $tilde(w)^((n+1))$ ist eine Dimension kleiner. Um diesen an
 das vorherige Ergebnis zu multiplizieren $w^((n+1)) = vec(0, tilde(w)^((n+1)))$
 sodass $Q^((n+1)) = mat(1, 0; 0, tilde(Q)^((n+1)))$ (alternativ: es wird auf
 $e^n$-Vektor gespiegelt)
 
 Insgesamt: in $M N^2-1/3 N^3$ Ops #h(2cm) $Q_N ... Q_2 Q_1 A =: R 
 #h(1cm)
 Q_1 ... Q_N =: Q$
 
 #grid(columns: 2*(1fr,), align: left,
 [
 === Lineares Ausgleichsproblem $norm(A x - b)_2 = min!$ 
 Berechne $A = Q R$ \
 Berechne $Q^T b = vec(c, d)$ \
 Löse $dash(R) x = c$ (Rückwärtssubstitution)
 ],
 [
 === Lösen von LGS $A x = b$ (falls $M=N$)
 $
 A x &= b \
 <=> Q R x &= b \
 <=> R x &= Q^T b \
 <=> R x &= c " (Rückwärtssubstitution)"
 
 $
 ],
 )
 
 #v(-1cm)
 === QR-Algorithmus (Eigenwertberechnung)
 
 1. setze $A_0 = A$ und $k = 0$
 2. zerlege $A_k = Q_k R_k$
 3. berchne $A_(k+1) = R_k Q_k$, erhöhe k um 1 und gehe zu 2 #h(2cm) 
    $O(N^3)$ Ops pro Iteration
 
 $k -> oo: A_k -> R = mat(lambda_1, *, *, *; , lambda_2, *, *; ,, dots.down, *;
 0,,,lambda_N), Q_k -> I_N$ #h(2cm) Konvergenz: $abs(lambda_N)/abs(lambda_(N-1))$
 
 Benutze *Hessenbergform* für weniger Ops pro Iteration $O(N^2)$. Kann durch N-2
 Householder-Transformation berechnet werden $Q^T A Q = H$
 
 Householder-Transformation $tilde(Q)_2 = I_(N-1) - 2 tilde(w)^2 (tilde(w)^2)^T
 in bb(R)^((N-1) times (N-1))$ mit $tilde(Q)_2 vec(a_(21), a_(31), dots.v, a_(N
 1)) = vec(*, 0, dots.v. 0) in bb(R)^(N-1), Q_2 = I_N - 2 w w^T =
   mat(1,;,tilde(Q)_2), w^2 = vec(0, tilde(w)^2)$
 
 #v(-.5cm)
 $
 A^((1)) = Q_2 A Q_2 \
 Q_(N-1) A^((N-3)) Q_(N-1) = Q_(N-1) ... Q_2 A Q_2 ... Q_(N-1) = Q^T A Q = H
 #h(1cm)
 , Q = Q_2 ... Q_(N-1)
 $
 
 Wenn wir eine Hessenbergmatrix QR-Zerlegung $H=Q R$, dann ist die Matrix
 $dash(H) = R Q$ ebenfalls Hessenbergmatrix, wegen $dash(H) = R Q = Q^T Q R Q =
 Q^T H Q$
 
 *QR-Algorithmus mit Shift* $mu_k = h_(NN)^((k))$ für schnellere Konvergenz
 $abs(tilde(h)_(n+1,n)^((k))) = O(abs((lambda_(n+1) - mu)/(lambda_n - mu)))$
 
 1. setze $H_0 = H, k = 0$ und wähle Toleranz $epsilon$
 2. zerlege $H_k - h_(N N)^((k')) I_N = Q_k R_K$
 3. berechne $H_(k+1) = R_k Q_k + h_(N N)^((k)) I_N$
 4. ist $abs(h_(N,N-1r^((k+1)))) <= epsilon (abs(h_(N N)^((k+1)) +
    abs(h_(N-1,N-1)^((k+1))))$ so akzeptiere $h_(N N)^((k+1))$ als EW. Sonst
    erhöhe k um 1 und gehe zu 2.
 
 
 #pagebreak()
 
 == Newton-Verfahren
 
 Problem: Für eine vorgegebene (nichtlineare) Funktion $f: D subset bb(R)^N ->
 bb(R)^N$ finde $x^* in D$ mit $f(x^*) = 0_N in bb(R)^N$
 
 Sei $x^0$ in der Nähe der Lösung $x^*$. Eine Taylor-Entwicklung liefert
 #v(-.8cm)
 $
    && 0 = f(x^*)                        = & f(x^0) + f'(x^0)(x^* - x^0) + overbrace(O(norm(x^0 - x^*)^2),"vernachlässigen" approx 0) \
 ~> &&0 =                                & f(x^0) + f'(x^0)(x^*-x^0) \
    &&" "=                               & f(x^0) + f'(x^0)x^* - f'(x^0)x^0 \
 <=>&& f'(x^0)x^0 - f(x^0)=              & f'(x^0) x^* \
 <=>&& f'(x^0)^(-1) (f'(x^0) x^0 - f(x^0)) = & x^* \
 <=>&& x^* approx                        & x^0 - f'(x^0)^(-1) f(x^0) \
  =>&& x^(k+1)                           = & x^k - f'(x^k)^(-1) f(x^k) \
    & ""#place()[Berechnung der inversen Jacobi-Matrix vermeiden] \
 <=>&& x^(k+1)                           = & -f'(x^k)^(-1) f(x^k) \           
 <=>&& f'(x^k) underbrace((x^(k+1) - x^k), =d^k) = & - f(x^k) \                      
 $
 
 1. Wähle Startwert $x^0 in D$ und Toleranz $epsilon$. Setze $k=0$
 2. Löse $f'(x^k)d^k = -f(x^k)$ (LGS durch LR-Zerlegung) \
    $x^(k+1) = x^k + d^k$
 3. Falls ($norm(d^k) < epsilon$): STOP \
    Sonst erhöhe k um 1 und gehe zu 2.
 
 Bem.: Man wählt nicht $norm(f(x^k)) < epsilon$ als STOP-Kriterium, weil das
 Problem $A f(x) = 0$ (A regulär) dann eine andere Lösung hätte. 
 $A f(x) = 0$ ändert Iterierten nicht, deshalb ist $norm(d^k) < epsilon$ besser,
 weil so bei $A f(x) = 0$ und $f(x) = 0$ die Lösung und Iterierten gleich bleiben. 
 
 ==== lokale quadratische Konvergenz
 
 Es gibt eine Kugel $kappa := kappa_p (x^*) = { x in bb(R)^N: norm(x-x^*) <= p }
 subset D$ sodass $x^*$ die einzige Nullstelle von f in $kappa$0ist. Zudem liegen
 die Folgeglieder $x^(k+1) = x^k - f'(x^;)^(-1) f(x^k)$ für jeden Startwert $x^0
 in kappa$ ebenfalls in $kappa$ und es gilt $lim_(k->oo) x^k = x^*$.
  Weiter existiert eine Konstante $c > 0$ mit $norm(x^* - x^k) <= c norm(x^* -
  x^(k-1))^2$ für $k in bb(N)$
 
 ==== Homotopie-Verfahren: guten Startwerd suchen
 $g: bb(R)^N -> bb(R)^N$ eine (einfache) Funktion mit bekannter Nullstelle
 $x^0$. Definiere $f_s (x) = (1-s) g(x) + s f (x)$, wählen $0 < s_1 < s_2 < ...
 < s_M = 1$ und lösen $f_(s_n) (x) = 0$ mit Newton mit Startvektor $x^(m-1)$ mit
 $f_(s_(m-1)) (x^(m-1) = 0$ für $m=1,...,M$
 
 ==== Dämpfungsstrategie: Verbesserung der Konvergenz
 
 Bestimme $t_k in (0,1]$ sodass für $x^(k+1) = x^k - t_k f'(x^k)^(-1) f(x^k)$ das
 Residuum kleiner wird $norm(f(x^(k+1))) < norm(f(x^k))$
 
 ==== Differenzenqoutient
 Wenn $f'(x^k)$ aufwendig zu berechnen, wähle $A_k approx f'(x^k)$ durch
 $(A_k)_(n m) = 1/h (f_n (x^k + h ^m) - f_n (x^k))$
 
 ==== Vereinfachtes Newton-Verfahren
 
 pro Iteration $f'$ auswerten und LR-Zerlegen ist teuer. Stattdessen $A approx
 f'(x^0)$ alse einmal $f'$ auswerten und LR-Zerlegen. $F(x) = x - A^(-1) f (x)$
 Nurnoch lineare Konvergenz. 
 
 ==== Globale Minimierung
 
 Wir suchen $x^* in bb(R)^N$ mit $f(x^*) <= f(x)$ für alle $x in bb(R)^N$. dh
 insbesodere grad $f(x^*) = 0_N$. Wir berechnen Nullstelle von $gradient f(x) =
 ("grad" f(x))^T$ mit Newton $x^(k+1) = x^k - H_f (x^k)^(-1) gradient f(x^k)$.
 $H_f$ Hessematrix. 
 
 ==== Nichlineare ausgleichrrechnung $norm(f(x))_2 = min!$
 
 Wähle Startwert $x^0 in bb(R)^N$ für $k=0,1,...$\
 1. löse lineares Ausglechsproblem $norm(f'(x^k) d^k + f(x^k))_2 = min$
 2. setze $x^(k+1) = x^k + d^k$
 
 #pagebreak()
 
 == Interpolation
 
 Problem: Suche zu den Stützstellen $(x_0,f_0), ..., (x_N, f_N)$ eine Funktion p
 mit $p(x_n) = f_n "für alle " n=0,...,N$
 
 === Polynom-Interpolation
 
 
 
 Problem: Suche zu N+1 Stützstellen $(x_0, f_0), ..., (x_0, f_0)$ ein Polynom p
 vom Grad $<= N (p in bb(P)_N)$, welches das Interpolationsproblem (s.o.) lösst:
 $p(x_n) = f_n " für alle " n = 0,...,N$, 
 
 ==== Eindeutigkeit des Interolationspolynoms zu einem Interpolationsproblem: 
 
 Nehmen wir an es gäbe zwei Polynome $p,q in bb(P)_N$ vom Grad $<= N$ welche das
 gleiche Interpolationsproblem lösen. $p(x_n) = q(x_n) = f_n "für " n = 0,...,N$.
 Dann wäre $p-q in bb(P)_N$ ein Polynom vom Grad $<= N$, welcehs $N+1$
 Nullstellen an den Stützstellen $x_n$ hat: $p(x_n) - q(x_n) = f_n - f_n = 0 "für alle " n=0,...,N$. 
 
 Nach dem Fundamentalsatz der Algebra kann ein Polynom vom Grad $=N$ nur maximal
 N Nullstellen haben außer es ist das konstante 0-Polynom $f=0$. 
 
 Damit das Polynom $p-q$ vom Grad $<= N $ also $N+1$ Nullstellen haben kann, muss
 nach dem Fundamentalsatz der Algebra $p - q = 0 <=> p = q$ das konstante
 0-Polynom sein, was bedeutet, das p und q identisch sein müssen. 
 
 ==== Lagrange-Basis
 
 Das n-te Lagrnge-Polynom $L_n$ zu den Stztzstellen $x_0, ..., x_N$ ist definiert
 als:
 
 $
 p(x) = sum_(n=0)^N f_n L_n (x) #h(1cm) 
 , L_n (x) = product_(m=0,m != n)^N (x-x_m)/(x_n-x_m) = sigma_(n m) = cases(1 ", falls "
 x_m = x_n, 0 ", falls " x_m != x_n) 
 
 $
 
 Lagrnge-Polynome $L_n$ bilden aufgrund von Skalarprodukt $angle.l p,q angle.r =
 sum_(n=0)^N p(x_n) q(x_n)$ eine orthonormalbasis und linearkombination
 ($bb(P)_N$).
 
 ===== Kondition
 
 Frage: Wie wirken sich Störungen der Eingabegrößen (Stützwerte $f_n$) auf die
 Lösung der Interpolation aus?
 
 $
 p(x) = sum_(n=0)^N f_n L_n (x) #h(2cm)
 tilde(p)(x) = sum_(n=0)^N tilde(f_n) L_n (x) \
 $
 $
 &abs(p(x) - tilde(p)(x)) = abs(sum_(n=0)^N (f_n - tilde(f_n)) L_n (x))
 <= sum_(n=0)^N abs(f_n - tilde(f_n)) abs(L_n (x)) 
 <= max_(n=0,...,N) abs(f_n - tilde(f_n)) sum_(n=0)^N abs(L_n (x)) \ 
 <=> &max_(x in [a,b]) abs(p(x) - tilde(p)(x)) <= underbrace( max_(x in [a,b]) sum_(n=0)^N
 abs(L_n (x)), =: Lambda_N >= 1) max_(n=0,...,N) abs(f_n - tilde(f_n))
 $
 
 Lebesgue-Konstante $Lambda_N$ ist invariant under Affinen Transformationen daher
 nur von relativer Lage der Stützstellen $x_n$ zueinander abhängig. 
 
 #let annotate(..args) = {
   box(place(bottom + center, dy: -10pt, stack(math.script(..args), line(angle: 90deg, length: 10pt))))
   sym.wj
   h(0pt, weak: true)
 }
 
 $
 #[==== Newton-Darstellung] #h(2cm)
 p(x) &= a_0 + a_1 w_1(x) + ... + #annotate([Leitkoeffizient])a_N w_N (x) #h(1cm)
 && w_n (x) = product_(m=0)^(n-1) (x - x_m)
 \
 "Aus Darstellung folgt "
 p(x) &= p_(0,N-1)(x) + a_N w_N (x) && p_(n,n)(x) = f_n = f_(n,n) \
 &= p_(0,N)(x) = f_(0,0) + f_(0,1) w_1(x) + ... + f_(0,N) w_N (x) #h(1cm) &&"mit " a_n = f_(0,n)
 $
 
 #grid(columns: 2, align: center + horizon,
 $
 #[*Lemma von Aitken*] #h(1cm)
 p_(n,k)(x) &= ((x_n - x) p_(n+1,k)(x) - (x_k-x)p_(n,k-1))/(x_n - x_k) \
 => "Für Koeffizienten" f_(n,k) &= (f_(,k-1) - f_(n+1,k)) / (x_n - x_k),
 #h(.5cm) "(Dividierte Differenzen)"
 $,
 $
 #h(.5cm) mat(delim: #none,
 f_(n,k-1);
 ,arrow.br;
 f_(n+1,k), arrow.r, f_(n,k);
 )
 mat(delim: #none,
 "2(N-1) Additionen";
 "N-1 Divisionen"
 )
 $
 )
 #grid(columns: 2,align: left+horizon, gutter: .5cm,
 [ *Interpolationsfehler* ],
 $
 f(x) - p(x) = w_(N+1)(x) (f^(N+1)(xi))/(N+1)! #h(1cm), 
 xi = xi(x) in (a,b) "Zwischenstelle" \
 $
 )
 #v(-.3cm)
 #grid(columns: 2,align: left+horizon, gutter: .5cm,
 [*Maximaler Interpolationsfehler/Approximationsgüte*],
 $
 max_(x in [a,b]) abs(f(x) - p(x)) <= max_(x in [a,b]) abs(w_(N+1)(x)) max_(x in [a,b]) abs(f^(N+1)(x))/(N+1)!
 $
 )
 
 ==== Tschebyscheff-Interpolation
 
 Tschebyscheff-Interpolationsformel
 
 $
 p(x) = 1/2 c_0 + c_1 T_1(x) + ... + c_N T_N (x)
 $
 
 mit Koeffizienten $c_m = 2/(N+1) sum_(n=0)^N f_n cos(m (2n +1)/(2N + 2) pi)$ mit
 $(N+1)^2$ Multiplikationen
 
 Tschebyscheff-Polynom vom Grad N ist auf $[-1,1]$ gegeben durch
 
 #grid(
 columns: 2, align: center + horizon, 
 $
 T_0(x) &= 1 \
 T_1(x) &= x \
 T_(N+1)(x) &= 2 dot T_N (x) - T_(N-1)(x) $,
 $
 "oder "
 T_N (x) = cos(N arccos(x))
 $
 )
 
 *Tschebyscheff-Stützstellen*: $x_n = cos((2n+1)/(2N + 2) pi) $ für alle $n=0,...,N$
 
 *Min-Max-Eigenschaft*: $max_(x in [-1,1]) abs(w_(N+1)(x))$ wird genau für
 Tschebyscheff-Stützstellen minimal mit $2^(-N)$
 
 
 mit FFT lässt sich Koeffizient $c_m$ in $O(N log N)$ Multiplikationen berechnen,
 wenn $N+1 = 2^k$ eine 2-er Potenz ist. 
 
 *Clenshaw-Algorithmus* zur Auswertung $p(x)$ mit N+2 Multiplikationen und 2N Additionen
 $
 d_(N+2) &= d_(N+1) = 0 \
 d_n &= c_n + 2x dot d_(n+1) - d_(n+2) " für " n=N,N-1,...,0 \
 "Dann gilt "& p(x) = 1/2 (d_0 - d_2)
 $
 
 #table(columns: 3,
 table.header([Zusammenfassung], [Bestimmung], [Auswertung]),
 [Monom],
 [LGS mit $1/3 N^3$ Ops],
 [Horner-Schema mit N Ops],
 [Lagrange],
 [direkt gegeben],
 [zu aufwendig],
 [Newton],
 [Dividierte Differenzen mit $N^2$ Ops],
 [Horner-Schema mit $2N$ Ops],
 [Tschebyscheff],
 [FFT mit $O(N log N)$ Ops],
 [Clenshaw-Algorithmus mit $2N$ Ops],
 )
 
 nur in Newton-Darstellung können zusätzliche Stützstellen hinzugefüngt werden,
 wobei nicht alle Koeffizienten nei berechnet werden müssen. 
 
 === Spline-Interpolation
 
 Problem: Suche zeimal stetig differenzierbare Funktion $s in C^2([a,b], bb(R))$
 mit mit $s(x_n) = y_n$ für $n=0,...,N$, soduss folgendes Integral minimal ist: 
 
 $
 integral_a^b abs(s''(x))^2 d x <= integral_a^b abs(s''(x) + epsilon h''(x))^2 d x
 \
 <=> integral_a^b s''(x)h''(x) d x 
 = sum_(n=1)^N integral_(x_(n-1))^(x_n) s''(x)h''(x) d x 
 = s''(b)h'(b) - s''(a)h'(a) =^! 0
 $
 für alle $epsilon in bb(R)$ und $C^2$-Funktionen $h: [a,b] -> bb(R)$ mit $h(x_n)
 = 0$ für alle $n=0,...,N$
 
 Aus dieser Gleichung ergeben sich derei Lösungen und damit drei Typen von Splines:
 
 #align(center, table(stroke: none, columns: 2, align: left+bottom,
 $s'(a) = v_0 " und " s'(b) = v_N$, [eingespannt / hermitisch],
 $ s''(a) = 0 = s''(b)$, [natürlich],
 $s'(a) = s'(b) " und " s''(a) = s''(b)$, [periodisch (sinnvoll für $overbrace(s(a) = s(b), y_0 = y_N)$)]
 ))
 
 not-a-knot-Bedingung falls Ableidung an Randpunkten unbekannt: $s'''_1(x_1) =
 s'''_2(x_1)$ und $s'''_(N-1) = s'''_N (x_(N-1))$
 
 #pagebreak()
 ==== kubischer interpolierender Spline
 
 Eine $C^2$-Funktion $s: [a,b] -> bb(R)$ heißt kubischer interpolierender Spline
 zu den Stützpunkten $(x_0, y_0), ..., (x_N, y_N)$ mit der Eigenschaft $a = x_0 <
 ... < x_N = b$, falls
 
 1. $s|_([x_(n-1),x_n]) in bb(P)_3$ für alle $n=1,...,N$
 2. $s(x_n) = y_n$ für alle $n=0,...,N$
 
 $
 s|_([x_(n-1),x_n]) = s_n (x) = 
 underbracket(y_(n-1) + (x - x_(n-1)) overbrace(y_(n-1,n), = (y_(n-1) - y_n)/(x_(n-1) - x_n)), "interpolierender Anteil") +
 underbracket((x-x_(n-1))(x-x_n)[alpha_n (x - x_(n-1)) + beta_n (x-x_n)],
 "glättender Anteil")
 $
 
 Daraus erbeben sich 2N Unbekannte $alpha_n, beta_n$. Diese erhält man durch die
 2N - 2 Glattheitsbedingungen + 2 Randbedingungen (s.o.). Also erhält man
 $alpha_n, beta_n$ aus LGS mit dim 2N. 
 
 $
 cases(reverse: #true, 
 s'_n (x_n) = s'_(n+1)(x_n),
 s''_n (x_n) = s''_(n+1)(x_n)
 ) " für " n=1,...,N-1 => 2N-2 " Bedingungen"
 $
 
 Die Berechnung der Unbekannten können wir vereinfachen, indem wir die $alpha_n,
 beta_n$ durch N+1 Unbekannte Momente $gamma_n$ darstellen:
 
 $
 cases(reverse: #true,
 s''_n (x_(n-1)) &= gamma_(n-1),
 s''_n (x_n) &= gamma_n
 ) " damit " s''_n (x_n) = s''_(n+1) (x_n) = gamma_n
 $
 
 Dies sind 2N Gleichungen für die N+1 Unbekannten $gamma_0, ..., gamma_N$. Durch
 den Ansatz sind N-1 Glattheitsbedingungen an die zweite Ableitung automatisch
 erfüllt. 
 
 Mit den Gitterweiten $h_n := x_n - x_(n-1)$ erhalten wir
 $#h(1cm) alpha_n = 1/(6 h_n) (gamma_(n-1) + 2 gamma_n) #h(1cm)
 beta_n = -1/(6 h_n) (2 gamma_(n-1) + gamma_n)$
 
 Eingesetzt in $s'$ in der Glattheitsbedingung $s'_n (x_n) = s'_(n+1) (x_n)$ folgt
 
 $
 h_n / 6 (gamma_(n-1) + 2 gamma_n) + h_(n+1) / 6 ( 2 gamma_n + gamma_(n+1)) =
 gamma_(n,n+1) - gamma_(n-1,n) =: d
 $
 
 mit Einschränkung auf eingespannte kubische Splines $s'_1(x_0) = 0_0, s'_N (x_N)
 = v_N$ erhalten wir außerdem
 
 $
 h_1/6 ( 2 gamma_0 + gamma_1) = y_(0,1) - v_0 =: d_0 
 #h(2cm)
 h_N / 6 (gamma_(N-1) + 2 gamma_N) = v_N - y_(N-1,N) =: d_N
 \
 "Insgesamt ein LGS: "
 1/6 underbrace( mat(
 2 h_1, h_1;
 h_1, 2(h_1 + h_2), h_2;
 ,dots.down, dots.down, dots.down;
 ,,h_(N-1), 2(h_(N-1) + h_N), h_N;
 ,,,h_N, 2h_N;
 ), #[
 strikt diagonaldominante symmetrische Tridiagonalmatrix (spd, regulär) \
 $=>$ mit Cholesky in $O(N)$ Ops lösbar. 
 ])
 vec(gamma_0, gamma_1, dots.v, gamma_(N-1), gamma_N)
 =
 vec(d_0, d_1, dots.v, d_(N-1), d_N)
 $
 
 ===== Kondition
 
 $tilde(s)(x)$ der Spline zu gestörten Daten $(x_n, tilde(y)_n)$, Lagrange-Spline
 $l_n (x_m) = cases(1 ", falls " n=m, 0 ", falls" n != m), #h(1cm) l'_n (a) = 0, l'_n (b) =0$
 
 $
 abs(s(x) - tilde(s)(x)) <= sum_(n=0)^N abs(y_n - tilde(y)_n) norm(l_n (x)) 
 <= sum_(n=0)^N abs(l_n (x)) max_(m=0,...,N) abs(y_m - tilde(y)_m) 
 <= Lambda_N max_(n=0,...,N) abs(x_m - tilde(y)_m)
 $
 
 $
 "mit " Lambda_N := max_(x in [a,b]) sum_(n=0)^N abs(l_n (x)) " also " 
 max_(x in [a,b]) abs(s(x) - tilde(s)(x)) <= Lambda_N max_(m=0,...,N) abs(y_m - tilde(y)_n) 
 $
 
 Bei äquidisanten Unterteilung gilt $Lambda_N <= 2$ für alle $N in bb(N)$
 
 #pagebreak()
 
 
 == Numerische Integration
 
 Problem: Berechne für eine gegebene Funktion $f: [a,b] -> bb(R)$ das
 Riemann-Integral. 
 
 === Eigenschaften von Integralen
 
 #table(columns: 2*(auto, 1fr,), stroke: none, 
 [Additiv:],[ $ integral_a^b f(x) d
 x = integral_a^c f(x) d x + integral_c^b f(x) d x$],
 [Linear:],[ $I(lambda f + mu g) = lambda I(f) + mu I(g)$],
 )
 Monoton: falls $f >= g$ auf $[a,b]$ dann $integral_a^b f(x) d x >= integral_a^b g(x) d x$ 
   \ $abs(I(f)) <= I(abs(f)) <=> cond_1(f) = I(abs(f)) / abs(I(f)) >= 1$
 
 Norm: $norm(f)_1 = I(abs(f))$
 
 === Quadraturformeln
 
 $
 integral_a^b f(x) d x 
 approx integral_a^b p(x) d x 
 = integral_a^b sum_(n=0)^N f(x_n) L_n (x) 
 = sum_(n=0)^N overbrace(integral_a^b L_n (x) d x, =(b-a)b_(n+1))
 f underbrace((x_n), #box()#place(horizon+center, $script(=a + c_(n+1)(b-a))$)) 
 = (b-a) sum_(k=1)^s b_k f( a + c_k (b-a))
 $
 
 #table(columns: 3,
 [Rechteckregel], 
 $s = 1, b_1 = 1, c_1 = 0, p = 1$,
 $I(f) approx (b - a) f(a)$,
 [Mittelpunktregel], 
 $s = 1, b_1 = 1, c_1 = 1/2, && p = 2$, 
 $I(f) approx (b-a) f((a+b)/2)$,
 [Trapezregel],
 $s = 2, b_1 = b_2 = 1/2, c_1 = 0, c_1 = 1, p = 2$, 
 $I(f) approx (b-a) (f(a) + f(b))/2$,
 [Simpsonregel],
 $s = 3, b_1 = b_3 = 1/6, b_2 = 4/6, c_1 = 0, c_2 = 1/2, c_3 = 1, p = 4$,
 $I(f) approx (b-a) 1/6 (f(a) + 4 f((a+b)/2) + f(b))$,
 )
 
 Eine QF $(b_k, c_k)_(k=1,...,s)$ besitzt die Ordnung p, falls sie für alle
 Polynome vom Grad $<= p - 1$ das Integral exakt berechnet, wobei p maximal ist.
 
 #TODO[BEWEIS ZU ORDNUNGBEDINGUNGEN!!!!!]\
 
 *Ordnungsbedingungen*: Eine QF besitzt genau dann die Ordnung p, falls 
 
 $
 1 / q = sum_(k=1)^s b_k c_k^(q-1) "für alle " q = 1, ..., p " aber nicht mehr für " q = p + 1 " gilt"
 $
 
 Für eine QF mit vorgegebenen Knoten $c_1 < ... c_s$ können die Gewichte genau
 dann mit $b_k = integral_0^1 underbrace(L_k (x), #place(center, $L_k (x) = product_(j=1,j != k)^s
 (x-c_j)/(c_k - c_j)$)) d x$ eindeutig bestimmt werden wenn $p >= s$. 
 
 *Symmetrische QF*: 
 
 $
 b_k = b_(s+1-k) #h(1.5cm) ,c_k = 1 - c_(s+1-k) "also symmetrisch um "1/2" verteilt"
 $
 
 Die Ordnung p einer symmetrischen QF ist gerade. #TODO[Beweis]
 
 Für symmetrische Knoten $c_k$ und $p >= s$ sind die eindeutigen Gewichte $b_k$
 ebenfalls symmetrisch.
 
 === QF mit erhöhter Ordnung
 
 Die maximale Ordnung einer QF ist 2s (insb. bei Gauß-QF)
 
 Bew.: $angle.l M,M angle.r = integral_0^1 M(x)^2 d x > 0, deg(M) = s$
 
 *Gauß-QF*: Ordnung 2s gegeben durch $c_k = 1/2(1+gamma_k)$ mit
 $gamma_1,...,gamma_s$ die Nullstellen des Legendre-Polynoms vom Grad s.
 Die Gewichte $b_k$ einer Gauß-QF sind positiv. Gauß-QF sind symmetrisch, weil
 Legendre-Nullstellen $gamma_1, ..., gamma_s$ symmetrisch zum Punkt 0 im Interval
 $(-1, 1)$ und so Knoten symmetrisch zu Punkt $1/2$ im Interval $(0,1)$. 
 
 *Summierte QF*: Um Fehler zu verkleinern zerlege $[a,b]$ in Teilintervalle
 $[x_(n-1), x_n], n=1,...,N$ mit Intervallängen $x_n - x_(n-1) = h_n$
 
 $
 integral_a^b f(x) d x approx  sum_(n=1)^N h_n sum_(k=1)^s b_k (x_(n-1) + c_k h_n)
 $