KIT-BOOKS/numerik/summary.typ

#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)_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)
$