diff --git a/numerik/summary.typ b/numerik/summary.typ new file mode 100644 index 0000000..b40e40e --- /dev/null +++ b/numerik/summary.typ @@ -0,0 +1,653 @@ +#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]) + +== 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)_1 = max_(n=1,...,N) sum_(m=1)^N abs(a_"nm") $, +$ 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_(NN)^((k')) I_N = Q_k R_K$ +3. berechne $H_(k+1) = R_k Q_k + h_(NN)^((k)) I_N$ +4. ist $abs(h_(N,N-1r^((k+1)))) <= epsilon (abs(h_(NN)^((k+1)) + + abs(h_(N-1,N-1)^((k+1))))$ so akzeptiere $h_(NN)^((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) = sum_(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 integral_a^b L_n (x) d x f(x) = (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 mi $b_k = integral_0^1 underbrace(L_k (x), #place(center, $L_k (x) = sum_(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) && 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) +$