\chapter{Numerische Integration}\label{KAPITEL-05}
%
Aufgabe der numerischen Integration (Quadratur) ist die Berechnung
von 
%
\[
 \int\limits^b_a f(x) dx, \quad a,b < \infty
\]
%
Typischerweise ist $f$ nicht direkt integrierbar.

\underline{Idee:} Approximiere $f(x)$ durch Funktion, die einfach zu
integrieren ist; Polynome, trigonometrische Polynome, 
Splines, rationale Funktionen, Exponentialsummen, \dots
%
\section{Newton--Cotes Formeln}
%
Die Idee der Newton--Cotes Formeln ist $f(x)$ durch ein Polynom zu interpolieren
und dieses zu integrieren. Bilde \"aquidistante Unterteilung von $[a,b] 
\quad x_i = a+ih  \quad i = 0,\cdots,n,$ mit $h = \frac{b-a}{n} \,\, n > 0 
\,\, n \in \N$ und bilde Interpolationspolynom in $\Pi_n$, das
%
\eq{5.1.1}
 P_n (x_i) = f(x_i) \quad i = 0, \cdots, n
\en
%
erf\"ullt. Satz (\ref{4.1.1}) $\Longrightarrow$
%
\eq{5.1.2}
 P_n (x) = \sum\limits^n_{i=0} f_i L_i (x).
\en
%
$\Longrightarrow$
%
\eqa\nonumber
 \int\limits^b_a P_n (x) dx & = & \sum\limits^n_{i=0} f_i \int\limits^b_a
 L_i (x) dx \quad \mbox { mit } x = a+hs\\
%
 & = & h \sum\limits^n_{i=0} f_i \int\limits^n_0 \Pi^n_{k=0 \atop k \neq i}
\label{5.1.3} \frac{s-k}{i-k} ds\\ 
%
 & = & h \sum\limits^n_{i=0} f_i \alpha_i = \frac{b-a}{n}
 \sum\limits^n_{i=0} \sigma_i f_i \quad \sigma_i \in \Z.\nonumber
\ena
%
Gewichte $\alpha_i$ sind unabh\"angig von $f_i, a,b$ jedoch abh\"angig von $n$.\\
Die Gewichte liegen tabelliert vor. Diese Formeln hei{\ss}en Newton--Cotes--Formeln.\\
Fehler:
%
\eq{5.1.4}
 \int\limits^b_a \left(P_n (x) - f(x) \right) dx = h^{p+1}
 \kappa \cdot f^{(p)} (\xi), \xi \in (a,b).
\en
%
\begin{tabelle}\label{5.1.5}

\begin{tabular}{c|l|r|l|l|}
 $n$ & $\sigma_i$ & $ns$ & Fehler & Name \\[0.3cm]\hline

 $1$ & $1 \,\, 1$ & $2$ & $\frac{h^3}{12} f^{(2)} (\xi)$ & Trapezregel \\[0.3cm]
%
 $2$ & $1 \,\, 4 \,\, 1$ & $6$ & $\frac{h^5}{90} f^{(4)} (\xi)$ &
 Simpsonregel\\[0.3cm]
%
 $3$ & $1 \,\, 3 \,\, 3 \,\, 1$ & $8$ & $\frac{h^5 3}{80} f^{(4)} (\xi)$ & 
 $3/8$ Regel \\[0.3cm]
%
 $4$ & $7 \,\, 32 \,\, 12 \,\, 32 \,\, 7$ & $90$ & $\frac{h^717 8}{945} f^{(6)} (\xi)$
 & Milne--Regel \\[0.3cm]
%
 $5$ & $19 \,\, 75 \,\, 50 \,\, 50 \,\, 75 \,\, 19$ & $288$ & $\frac{h^7 275}{12096}
 f^{(6)} (\xi)$ & -- \\[0.3cm]
%
 $6$ & $41 \,\, 216 \,\, 27 \,\, 272 \,\, 27 \,\, 216 \,\, 41$ & $840$ &
 $\frac{h^9 9}{1400} f^{(8)} (\xi)$ & Weddle--Regel
\end{tabular}
\end{tabelle}
%
gr\"o{\ss}ere $n \Longrightarrow $ negative Gewichte.\\
Wegen der schlechten Approximationseigenschaften (global) von Polynomen
werden die Regeln st\"uckweise auf Teilintervalle angewendet und dann aufsummiert.
%
\begin{example}[Summierte Trapezregel: st\"uckweise lineare 
Interpolation]\label{5.1.6} 
\begin{center}
\unitlength 1cm
\begin{picture}(8,8)
% \put(0,0){\epsfig{file=includeps/trapez.ps,width=8cm,clip=}}
\put(0,0){\includegraphics[width=0.60\linewidth,height=0.50\linewidth]{includeps/trapez}}
\end{picture}
\end{center}
\end{example}
 
%
Teilintervalle $[x_i, x_{i+1}], x_i = a+ih \quad i = 0, \cdots, N \quad
h = \frac{b-a}{N}.$ Dann ist die summierte Trapezregel
%
\eqa\label{5.1.7}
 T(h) & := & \sum\limits^{N-1}_{i=0} \frac{h}{2} \left[ f(x_i) + f(x_{i+1})
 \right] \\
%
 & = & h \left[ \frac{f(a)}{2} + f(a+h) + \cdots + f(b-h) +
 \frac{f(b)}{2} \right] \nonumber
\ena
%
F\"ur jedes Teilintervall ist der Fehler
%
\eq{5.1.8}
 \frac{h}{2} \left[ f(x_i) - f(x_{i+1}) \right] - 
 \int\limits^{x_{i+1}}_{x_i} f(x) dx = \frac{h^3}{12} f^{(2)} (\xi_i)
 \mbox { f\"ur } \xi_i \in (x_i, x_{i+1}).
\en
%
Insgesamt ergibt sich also 
%
\eqa\label{5.1.9}
 \left| T(h) - \int\limits^b_a f(x) dx \right| & = & \left| \sum\limits^{n-1}_{i=0} \frac{h^3}{12}
 f^{(2)} (\xi_i) \right| \\
 & \le & \frac{h^3}{12} N \cdot \underbrace{\sup\limits_{\xi \in [a,b]}f''(\xi)}_{=M}
  = \frac{b-a}{12} h^2 M. \nonumber
\ena
%
Fehler geht mit $h^2$ gegen $0$, \underline{Verfahren 2.~Ordnung}.
%
\begin{example}[Mittelpunktregel]\label{5.1.10}~
\begin{center}
\unitlength 1cm
\begin{picture}(8,8)
% \put(0,0){\epsfig{file=includeps/mittelp.ps,width=8cm,clip=}}
\put(0,0){\includegraphics[width=0.60\linewidth,height=0.50\linewidth]{includeps/mittelp}}
\end{picture}
\end{center}
\end{example}

%
\begin{example}\label{5.1.11}
Sei $N$ gerade und verwende Simpsonregel auf $[x_{2i}, x_{2i+1}, x_{2i+2}],\;
i = 0,\cdots,\frac{N}{2} -1.$ Der N\"aherungswert auf jedem der Teilintervalle
ist 
\[
 \frac{h}{3} \left[ f(x_{2i})+4f (x_{2i+1})+f(x_{2i+2}) \right].
\]
%
Summation ergibt
%
\eqa\label{5.1.12}
 S(h) & = & \frac{h}{3} \left[ (f(a) + 4f (a+b) + 2f(a+2h)+4f(a+3h) + \right. \\
  & & \left. + 2f (b-2h) + 4f (b-h) + f(b) \right] \nonumber
\ena
%
und als Fehler
%
\eqa\label{5.1.12'}
 \left| S(h) - \int\limits^b_a f(x) dx \right| & \le & 
 \frac{h^5}{90} \cdot \frac{N}{2} \cdot
 \underbrace{\sup\limits_{\xi \in [a,b]} f^{(4)} (\xi)}_{=M} \\
& & = \frac{b-a}{180} h^4 \cdot M. \nonumber
\ena
\underline{Verfahren 4.~Ordnung}.
\end{example}
%
Weitere Quadraturformeln erh\"alt man durch andere Approximationsaufgaben
f\"ur $f(x)$.\\
Zum Beispiel \underline{Gauss--Quadratur} f\"ur Integrale der Form
\[
 \int\limits^b_a w(x) f(x) dx \approx \sum\limits^n_{i=1} w_i f(x_i).
\]
%
wobei $w_i, f_i$ so gew\"ahlt werden, da{\ss} der Fehler
\[
 \int\limits^b_a w(x) f(x) dx - \sum\limits^n_{i=1} w_i f(x).
\]
%
f\"ur Polynome m\"oglichst hohen Grades verschwindet $\longrightarrow$ \"Ubung.\\
Andere Ideen sp\"ater bei der L\"osung von Dgl.
%
\section{Extrapolation}
Untersuche das asymptotische Verhalten der Trapezsumme $T(h)$ f\"ur $h \to 0.$ Dazu
brauchen wir zuerst einige Eigenschaften der \underline{Bernoulli--Polynome}, die wie
folgt rekursiv definiert sind.
%
\begin{definition}\label{5.2.1}
Die Polynome $B_k(x), k = 0,1,2,\cdots,$ die durch die rekursive Folge
%
\eq{5.2.2.}
\begin{array}{l}
 B_0 (x) =1, B'_k (x) = k \cdot B_{k-1} (x) \quad k \ge 1,\\
 \int\limits^1_0 B_k (x) \, dx = 0,\; k \ge 1 \nonumber
\end{array}
\en
%
definiert sind, hei{\ss}en \underline{Bernoulli--Polynome.}\\
Eigenschaften der Bernoulli--Polynome:
%
\eq{5.2.3}
 B_k (x) = a_k + k \int\limits^x_0 B_{k-1} (t) \, dt \quad k \ge 1
\en
%
wobei $a_k$ so zu bestimmen ist, da{\ss} $\int\limits^1_0 B_k (x) \, dx = 0.$
Man erh\"alt sofort
%
\eq{5.2.4}
\begin{array}{l}
 B_1 (x) = x - \frac{1}{2}, B_2 = x^2 - x + \frac{1}{6}, 
 B_3(x) = x^3 - \frac{3}{2}x^2 + \frac{1}{2} x\\
 B_4 (x) = x^4 - 2x^3 + x^2 - \frac{1}{30}, \cdots \nonumber
\end{array}
\en
%
und es gilt das folgende Lemma.
\end{definition}
%
\begin{lemma}\label{5.2.5}~ 
\begin{itemize}
\item [i)] $B_k(0) = B_k(1)$ f\"ur $k \ge 2$.
\item [ii)] Die Polynome $P_k (x) = B_k (\frac{1}{2} + x)$ sind f\"ur
gerades $k$ gerade, d.h. $P(x) = P(-x)$ und f\"ur ungerades $k$ ungerade
Polynome, d.h. $P(x) = -P(-x)$ in $x$.
\item [iii)] $B_{2k+1}(0)=0$ f\"ur alle $k \ge 1.$
\end{itemize}
\end{lemma}
%
\proof
\begin{itemize}
\item [i)] Aus (\ref{5.2.3})
\[
 B_k (1) - B_k (0) = k \int\limits^1_0 B_{k-1} (t) \, dt = 0 
 \mbox { f\"ur } k \ge 2.
\]
%
\item [ii)] $P_0 (x) \equiv 1$ ist gerade. Sei $P_{2k}(x)$ f\"ur
$k \ge 0$ gerades Polynom. Zeige, dass $P_{2k+2}(x)$ wieder gerade ist.
Es gilt
%
\[
 P'_{2k+1} (x) = B'_{2k+1} (\frac{1}{2} + x) = (2k+1) B_{2k}
(\frac{1}{2}+x) = (2k+1) P_{2k} (x) 
\]
%
Also
\[
 \int\limits^{\frac{1}{2}}_{-\frac{1}{2}} P_{2k+1}'(x) dx = 0.
\]
Aus $P_{2k}(x)$ gerade folgt, da{\ss}
\[
 q(x) := (2k+1) \int\limits^x_0 P_{2k} (t) \, dt
\]
%
ein ungerades Polynom ist und es gilt
%
\[
 \int\limits^{\frac{1}{2}}_{-\frac{1}{2}} q(x) dx 
 = - \int\limits^{-\frac{1}{2}}_{+\frac{1}{2}} q (-x)\, dx =
 \int\limits^{-\frac{1}{2}}_{\frac{1}{2}} q(x) dx = 
 -\int\limits^{\frac{1}{2}}_{-\frac{1}{2}} q(x) dx = 0.
\]
%
Weiterhin ist 
$q'(x) = (2k+1) P_{2k}(x)$ und $\int\limits^{\frac{1}{2}}_{-\frac{1}{2}}
q(x) dx = 0$
$\Longrightarrow q (x) = P_{2k+1}(x)$ ungerade.\\
$P_{2k+2}(x) = (2k+2) \int\limits^x_0 P_{2k+1} (t) \, dt + a$ 
$\Longrightarrow P_{2k+2}$ ist wieder gerade. $\Longrightarrow$ ii).
%
\item[iii)] F\"ur $k \ge 1$ gilt
%
\[
\begin{array}{l}
 B_{2k+1} (1) = B_{2k+1}(0) = P_{2k+1}(-\frac{1}{2}) = - P_{2k+1} (\frac{1}{2})
 = - B_{2k+1}(1)\\
 \Longrightarrow B_{2k+1} (0) = B_{2k+1}(1) = 0.
\end{array}
\]
\end{itemize}
\eproof
%
$\beta_k = (-1)^{k+1} B_{2k} (0)$ hei{\ss}en \underline{Bernoullische Zahlen}.\\
$\beta_1 = \frac{1}{6}, \beta_2=\frac{1}{30}, \beta_3=\frac{1}{42}, 
\beta_4=\frac{1}{30}.$\\
Wir definieren nun neue Funktionen
%
\eq{5.2.6}
 S_k (x) := B_k(x-i) \mbox { f\"ur } x \in E_i := \{ x|i < x <i+1\}, \; i \le k.
\en
%
Diese Funktionen sind periodisch mit Periode 1 und durch die folgende
Rekursion bestimmt.
%
\eq{5.2.7}
 S_k (x) = S_k (0) + k \int\limits^x_0 S_{k-1} (t)\,dt,\; k \ge 2
\en
%
$S_1$ ist st\"uckweise linear mit Sprungstellen bei
$i \in \N$ und Nullstellen bei $i + \frac{1}{2}$ und es gilt wegen Lemma
\ref{5.2.5}
%
\eq{5.2.8}
\begin{array}{l}
 B_k (0) = S_k (0) = S_k (1) = S_k (2) = \cdots \\
 S_{2k+1} (0) = 0, S_{2k} (0) = (-1)^{k+1} \beta_k. \nonumber
\end{array} 
\en
%
Damit k\"onnen wir nun die Euler--MacLaurin'sche Summenformel beweisen.
%
\begin{theorem}\label{5.2.9}
 F\"ur $f \in {\cal C}^{2m+2} [a,b]$ hat $T(h)$ die Entwicklung
%
\eq{5.2.10}
 T(h) = \tau_0 + \tau_1 h^2 + \tau_2 h^4 + \cdots
 + \tau_m h^{2m} + \alpha_{m+1} (h) h^{2m+2}
\en
%
mit $\tau_0 := \int\limits^b_a f(x) \, dx.$ Dabei sind die $\tau_i$ von 
$h$ unabh\"angige Konstanten und $|\alpha_{m+1} (h) | \le M$
f\"ur alle $h = \frac{b-a}{n}, n \in \N.$
\end{theorem}
%
\proof
Betrachte $g(t) := f(a+th).$ Damit k\"onnen wir nun das Integral
\[
 \int\limits^n_0 S_1 (t) g' (t) \, dt
\]
%
durch partielle Integration umformen. $S_1(x)$ st\"uckweise linear
aus $B_1 (x) = x - \frac{1}{2}.$
\[
 \int\limits^1_0 S_1 (t) g' (t) \, dt = \int\limits^1_0 B_1 (t) g' (t) \, dt =
 \left.B_1 (t) g(t) \right|^1_0
 - \int\limits^1_0 B_0 (t) g(t) \, dt = \frac{1}{2} \left[ g(1) + g(0) \right]
 - \int\limits^1_0 g(t) \, dt.
\]
%
Analog
%
\begin{eqnarray*}
 \int\limits^{i+1}_i S_1 (t) g' (t) \, dt &=& \int\limits^1_0 B_1 (t) g' (t+i) \, dt\\
 &=& \frac{1}{2} \left[ g(i+1) + g(i) \right] - \int\limits^{i+1}_i g(t) \, dt,\;
 i = 0, \dots, n-1
\end{eqnarray*}
%
\eq{5.2.11}
 \Longrightarrow
  \frac{g(0)}{2} + g(1) + \cdots + g(n-1) + \frac{g(n)}{2} - 
  \int\limits_0^n g(t) \, dt
  = \int\limits^n_0 S_1 (t) g'(t) \, dt 
\en
%
F\"ur die rechte Seite gilt
%
\begin{eqnarray*}
  \int\limits^n_0 S_1 (t) g'(t) \, dt &=& \frac{1}{2!}\left( S_2 (t) g'(t) |^n_0 -
  \int\limits^n_0 S_2 (t) g'' (t) \, dt\right)\\
%
 &=& \frac{\beta_1}{2!} \left[ g' (n) - g'(0) \right]
   - \frac{1}{2!} \int\limits^n_0 S_2 (t) g'' (t) \, dt \\
%
 &=& \frac{\beta_1}{2!} \left[ g'(n)-g'(0) \right] - \frac{1}{3!}
   S_3 (t) g'' (t) |^n_0
 + \frac{1}{3!} \int\limits^n_0 S_3 (t) g^{(3)} (t) \, dt\\
& =& \frac{\beta_1}{2!} \left[ g'(n) - g'(0) \right] + \frac{1}{3!} 
   \int\limits^n_0 S_3 (t) g^{(3)}\, dt 
\end{eqnarray*}
%
$2m$ mal partielle Integration ergibt
%
\eq{5.2.12}
 \begin{array}{l}
 \frac{g(0)}{2} + g(1) + \cdots + g(n-1) + \frac{g(n)}{2} -
 \int\limits^n_0 g(t) \, dt\\
%
 = \sum\limits^m_{k=1} (-1)^{k+1} \frac{\beta_k}{(2k!)} 
   \left[ g^{(2k-1)} (n)- g^{(2k-1)} (0) \right] + R_{m+1} \nonumber
\end{array}
\en
%
mit 
%
\eqa\nonumber
 R_{m+1}& = &\frac{1}{(2m+1)!} \int\limits^n_0 S_{2m+1} (t) g^{(2m+1)} (t) \, dt\\
\nonumber  & =&~ \frac{1}{(2m+2)!}\left.(S_{2m+2} (t)-S_{2m+2}(0))g^{(2m+1)}(t)\right|_0^n\\
\nonumber&&       - \frac{1}{(2m+2)!} \int\limits^n_0  
       \left[S_{2m+2} (t)-S_{2m+2}(0)   \right] g^{(2m+2)} (t) \, dt\\
\label{5.2.13}  & =& - \frac{1}{(2m+2)!} \int\limits^n_0  
       \left[S_{2m+2} (t)-S_{2m+2}(0)   \right] g^{(2m+2)} (t) \, dt
\ena
%
%
Es ist $\int\limits^n_0 g(t) \, dt = \frac{1}{h} \int\limits^b_a f(x) dx$ und
$g^{(k)} (t) = h^k f^{(k)} (a+th) \quad k = 0,1,2, \cdots$
\begin{eqnarray*}
\Longrightarrow 
\frac{g(0)}{2} + g(1) + \cdots + g(n-1) + \frac{g(n)}{2} &=&
\frac{f(a)}{2} + f(a+h) + \cdots + f(b-h) + \frac{f(b)}{2}\\ &=&
 \frac{1}{h} T(h).
\end{eqnarray*}
%
Also folgt
%
\eq{5.2.14}
 T(h) = \tau_0 + \tau_1 h^2 + \cdots + \tau_m h^{2m} + \alpha_{m+1} (h) h^{2m+2}
\en
%
mit $\tau_0 = \int\limits^b_a f(x) dx$.
%
\eq{5.2.15}
 \tau_k = \frac{(-1)^{k+1}}{(2k)!} \beta_k \left[ f^{(2k-1)} (b) - f^{(2k-1)} (a) \right]
\en
%
\eq{5.2.16}
 \alpha_{m+1} = \frac{-1}{(2m+2)!} \int\limits^b_a f^{(2m+2)} (x)
 \left[ S_{2m+2} \left( \frac{x-a}{n} \right) - S_{2m+2}(0) \right] \, dx
\en
%
Da $S_{2m+2}$ stetig und periodisch, folgt, da{\ss} es eine von $h$
unabh\"angige Schranke f\"ur $|\alpha_{m+1}^{(n)}|$ gibt.
\eproof
%




Die Euler--MacLaurinformel l\"a{\ss}t sich nun zur
\underline{Extrapolation} verwenden (Romberg/Bulisch).

Wenn man das Restglied vernachl\"assigt, ist $T(h)$ ein Polynom  in $h^2$,
das f\"ur $h=0$ den Wert $\tau_0 = \int\limits^b_a f(x)\, dx$ liefert.

Idee: Bestimme f\"ur verschiedene Schrittweiten\\
$h_0 = b-a, h_1 = \frac{h_0}{n_1}, \cdots, h_m = \frac{h_0}{n_m}, \; n_i \in \N,$
bestimme die zugeh\"orige Trapezsumme
%
\eq{5.2.17}
 T_{i,0} = T(h_i),\; i = 0, \cdots, m
\en
%
und dann durch Interpolation das Polynom
%
\eq{5.2.18}
 \tilde{T}_{m,m} (h) = a_0 + a_1 h^2 + \cdots + a_m h^{2m},
\en
%
da{\ss} die Interpolationsbedingungen
%
\eq{5.2.19}
 \tilde{T}_{m,m} (h_i) = T(h_i), \; i=0,\dots,m
\en
%
erf\"ullt. Werte dann dies Polynom bei $h=0$ mit dem Neville--Schema aus, dies
liefert i.a. sehr gute N\"aherung f\"ur Integral.\\
Sei $\tilde{T}_{i,k} (h)$ das Polynom vom Grad $k$ in $h^2$, da{\ss}
%
\eq{5.2.20}
 \tilde{T}_{i,k} (h_j) = T_{j,0}, \; j = i-k, \cdots,i
\en
%
erf\"ullt. Damit gilt f\"ur die extrapolieren Werte $T_{i,k}:= \tilde{T}_{i,k} (0).$
%
\eq{5.2.21}
 T_{i,k} = T_{i,k-1} + \frac{T_{i,k-1} - T_{i-1,k-1}}{\left(\frac{h_{i-k}}{h_i}\right)^2-1},\;
 1 \le k \le i \le m
\en
%
$\longrightarrow$ \"Ubung \ref{numI/uebung069}
%
$\left\{\begin{array}{l}
  h_0 = b-a, h_1 = \frac{h_0}{2}, \cdots \Longrightarrow T_{11}
 \mbox { Simpsonregel}. \\
  h_0 = b-a, h_1 = \frac{h_0}{3}, \cdots \Longrightarrow T_{11} ~ 
 \frac{3}{8}-\mbox { Regel}. 
\end{array}\right.
$
%

Extrapolationsfolgen
\begin{itemize}
\item[a)] $h_0 = b-a, h_1 = \frac{h_0}{2}, h_2 = \frac{h_1}{2}, h_3= \frac{h_2}{2},
\cdots$ (Romberg Folge)
\item[b)] $h_0 = b-a, h_1 = \frac{h_0}{2}, h_2 = \frac{h_0}{3}, h_3 = \frac{h_1}{2},
h_4 = \frac{h_2}{2}, h_5 = \frac{h_3}{2}, \cdots$ (Bulirschfolge).
\end{itemize}
b) besser, da Rechenaufwand weniger schnell ansteigt. Vorherige Werte k\"onnen
verwendet werden.
%
$\longrightarrow$ \"Ubung \ref{numI/uebung069}.

\underline{Fehlerbetrachtung:} F\"ur ein $\xi \in [a,b]$
%
\eq{5.2.22}
 T_{mm} - \int\limits^b_a f(x) \, dx = -\frac{1}{(2m+2)!} f^{(2m+2)} (\xi)
 \int\limits^b_a K(x) \, dx
\en
%
mit
%
\eq{5.2.23}
 K(x) = \sum\limits^m_{i=0} c_{m_i}h_i^{2m+2} \left[ S_{2m+2}
 \left(\frac{x-a}{h_i}\right)  - S_{2m+2} (0) \right]
\en
%
%
\eq{5.2.24}
\Longrightarrow
 T_{mm} - \int\limits^b_a f(x) dx = (b-a) h_0^2 \cdots h_m^2
 \frac{\beta_{m+1}}{(2m+1)!}
 f^{(2m+2)} (\xi)
\en
%
Details siehe z.B. Stoer, {\it Numerische Mathematik}.




  
\begin{small}
\input{kap5u}
\end{small}