test

Wir wollen nun Integrale berechnen. Oft lässt sich ${\int }_a^{b}f(x)dx$ nicht analytisch lösen (z.B. elliptisches Integral beim Pendel) oder wir kennen $f$ nur an einzelnen Punkten → numerische Approximation.

Note: Wir stellen uns vor, dass jede Auswertung von $f$ teuer ist → wir wollen mit möglichst wenigen Funktionsauswertungen möglichst genau sein.
In Python: scipy.integrate.quad (adaptiv, QUADPACK). Selbst implementiert werden meist Trapez- und Simpson-Regel.

5.3 Grundbegriffe und -ideen

5.3.1 Quadraturformel

Quadraturformel (Def. 5.3.1)

Das Integral wird durch eine gewichtete Summe von Funktionswerten approximiert:

$$I(f):={\int }_a^{b}f(x)dx\approx \sum _{i=1}^n{w}_{i}f(c_i)=:Q_n(f;a,b)$$

wobei $w_i$ die Gewichte und $c_i\in [a,b]$ die Knoten der Quadraturformel sind.

Fehler (Def. 5.3.2)

$$E(n):=\left|{\int }_a^{b}f(x)dx-Q_n(f;a,b)\right|$$

Algebraische Konvergenz: $E(n)=O(n^{-p})$, $p>0$. Exponentielle Konvergenz: $E(n)=O(q^n)$, $0\le q<1$.

Note: Das ist die Riemann-Summe aus der Analysis, nur dass wir jetzt frei wählen, wo wir auswerten ($c_i$) und wie stark wir jeden Wert gewichten ($w_i$).

• links, rechts, Mitte, … → das sind verschiedene Wahlen von $c_i$
• die Wahl der $c_i$ und $w_i$ entscheidet, wie gut die Formel ist

5.3.2 Idee: Ersetze $f$ durch etwas Einfacheres

Wir approximieren $f$ durch eine Funktion, die wir exakt integrieren können → meistens ein Polynom $p$ mit $p(c_i)=f(c_i)$ (Interpolation, Kap. 2):

$${\int }_a^{b}f(x)dx\approx {\int }_a^{b}p(x)dx$$

Repetition Lagrange-Basis (Kap. 2.3): Zu Knoten $x_0,\dots ,x_n$ sind die Lagrange-Polynome

$${\ell }_i(x)=\prod _{j\ne i}\frac{x-x_j}{x_i-x_j},{\ell }_i(x_j)=\{\begin{array}{ll}1,& j=i\\ 0,& j\ne i\end{array}$$

Sie bilden eine Basis von den Polynomen vom Grad $\le n$, und das Interpolationspolynom ist einfach

$$p(x)=\sum _{j=0}^{n}f(x_j){\ell }_j(x)$$

Note: Kardinalbasis = eine Basis mit genau dieser Eigenschaft ${\ell }_i(x_j)={\delta }_{ij}$ (”$1$ an meinem Knoten, $0$ an allen anderen”).
→ die Koeffizienten in dieser Basis sind direkt die Funktionswerte $f(x_j)$, man muss kein Gleichungssystem lösen (anders als bei der Monombasis / Vandermonde-Matrix).

• wie $e_1,\dots ,e_n$ im ${\mathbb{R}}^n$: der $j$-te Eintrag eines Vektors ist seine Koordinate.

Quadraturgewichte aus den Lagrange-Polynomen (Bem. 5.3.7)

Integriert man das Interpolationspolynom, erhält man eine Quadraturformel mit Knoten $x_j$ und Gewichten

$$w_j={\int }_a^b{\ell }_j(x)dx,j=0,\dots ,n$$

Diese Formel ist exakt für alle Polynome vom Grad $\le n$.

Proof: Linearität des Integrals:

$${\int }_a^{b}p(x)dx={\int }_a^b\sum _{j}f(x_j){\ell }_j(x)dx=\sum _{j}f(x_j)\underset{=:w_j}{\underbrace{{\int }_a^b{\ell }_j(x)dx}}$$

Ist $f$ selbst ein Polynom vom Grad $\le n$, dann ist $p=f$ (Eindeutigkeit der Interpolation) → exakt.

Note: Warum ist das Lagrange? / Sind die $w_i$ durch die Knoten festgelegt? → Ja. Fixiert man $n+1$ Knoten und verlangt, dass die Formel alle Polynome vom Grad $\le n$ exakt integriert, dann gibt es genau eine Wahl der Gewichte, und zwar $w_j=\int {\ell }_j$:

• setze $f={\ell }_k$ ein: $Q({\ell }_k)={\sum }_j{w}_j{\ell }_k(x_j)=w_k$ (Kardinaleigenschaft!) und das muss $=\int {\ell }_k$ sein.
• die Gewichte sind also nicht die “Polynomkoeffizienten”, sondern die Integrale der Basisfunktionen. Die “Koeffizienten” sind die $f(x_j)$.

Fehlerschranke (Bem. 5.4.1)

Für $f\in C^n([a,b])$ und Interpolation mit $n$ Knoten (Grad $n-1$):

$$\left|{\int }_a^{b}f(t)dt-Q_n(f)\right|\le \frac{1}{n!}(b-a{)}^{n+1}\Vert f^{(n)}{\Vert }_{L^{\infty }([a,b])}$$

→ Zwei Dinge sind wichtig: die Glattheit von $f$ und die Länge $b-a$ des Intervalls. Kurzes Intervall = kleiner Fehler → nächste Idee.

5.3.3 Teile und Herrsche (summierte Quadratur)

Wir teilen $[a,b]$ in $N$ Teilintervalle der Länge $h=\frac{b-a}{N}$, $x_k=a+kh$, und wenden auf jedem eine einfache Formel an:

$${\int }_a^{b}f(x)dx=\sum _{k=0}^{N-1}{\int }_{x_k}^{x_{k+1}}f(x)dx\approx \sum _{k=0}^{N-1}Q(f;x_k,x_{k+1})$$

Note: Lokaler Fehler $O(h^{n+1})$, davon gibt es $N=\frac{b-a}{h}$ Stück → globaler Fehler $O(h^n)$. Beim Aufsummieren geht ein Faktor $h$ verloren.

Note: Man muss die Teilintervalle nicht gleich gross machen: dort, wo $f$ stark variiert (gross $|f^{\prime \prime }|$), feiner unterteilen → adaptive Quadratur (Kap. 5.6).

5.3.4 Referenzintervall

Man definiert Quadraturformeln einmal auf $[-1,1]$ (Knoten ${\hat{c}}_j$, Gewichte ${\hat{w}}_j$) und bildet sie affin auf jedes $[a,b]$ ab:

$${\int }_a^{b}f(t)dt=\frac{1}{2}(b-a){\int }_{-1}^1\hat{f}(\tau )d\tau ,\hat{f}(\tau ):=f\left(\frac{1}{2}(1-\tau )a+\frac{1}{2}(1+\tau )b\right)$$ $$c_j=\frac{1}{2}(1-{\hat{c}}_j)a+\frac{1}{2}(1+{\hat{c}}_j)b,w_j=\frac{1}{2}(b-a){\hat{w}}_j$$

5.3.5 Ordnung

Ordnung (Def. 5.3.8)

Eine Quadraturformel hat Ordnung $n+1$, wenn sie alle Polynome vom Grad $\le n$ exakt integriert (und $x^{n+1}$ nicht mehr).

Note: Ordnung $d$ $\iff$ exakt bis Grad $d-1$. Hat eine Formel Ordnung $2$, ist sie auch für Grad $0$ exakt; man gibt immer das grösste solche $d$ an.

• Warum "$+1$"? Ein Polynom vom Grad $d-1$ hat $d$ Koeffizienten → Ordnung = Anzahl Freiheitsgrade, die man exakt trifft.
• Prüfen in der Praxis: setze $1,x,x^2,\dots$ ein, das erste falsche Monom $x^d$ gibt Ordnung $d$.

5.3.6 Symmetrische Formeln

Symmetrisch (Def. 5.3.9)

Eine Quadraturformel auf $[-1,1]$ heisst symmetrisch, falls

$$w_i=w_{n+1-i}\text{und}{c}_i=-c_{n+1-i}$$

(Knoten spiegelbildlich um $0$, gespiegelte Knoten haben gleiche Gewichte).

Ordnung symmetrischer Formeln (Satz 5.3.11)

Die Ordnung einer symmetrischen Quadraturformel ist gerade.

Proof (Idee): Angenommen, die Formel ist exakt bis Grad $2m-2$. Sei $f\in {\mathcal{P}}_{2m-1}$, schreibe $f(x)=a{x}^{2m-1}+g(x)$ mit $\mathrm{deg}g\le 2m-2$.

• $g$ wird nach Annahme exakt integriert.
• $x^{2m-1}$ ist ungerade: ${\int }_{-1}^1{x}^{2m-1}dx=0$.
• Die Formel liefert auch $0$: Paare $(c_i,-c_i)$ mit gleichem Gewicht geben $w_i{c}_i^{2m-1}+w_i(-c_i{)}^{2m-1}=0$ (und ein Knoten bei $0$ gibt $0$).
• Linearität von $Q$ → $Q(f)=a\cdot 0+Q(g)=\int f$.

Also bekommt man Grad $2m-1$ gratis → die Ordnung springt von $2m-1$ auf $2m$, ist also gerade.

5.4 Äquidistante Stützstellen

Die drei wichtigsten Formeln (alle symmetrisch, Bem. 5.3.10), mit $m=\frac{a+b}{2}$:

	Mittelpunkt	Trapez	Simpson
Knoten	$m$	$a,b$	$a,m,b$
Gewichte	$b-a$	$\frac{b-a}{2},\frac{b-a}{2}$	$\frac{b-a}{6},\frac{4(b-a)}{6},\frac{b-a}{6}$
interpoliert durch	Konstante	Gerade	Parabel
exakt bis Grad	$1$	$1$	$3$
Ordnung	$2$	$2$	$4$
summiert: Fehler	$O(h^2)$	$O(h^2)$	$O(h^4)$

5.4.1 Mittelpunkt-Regel

Mittelpunkt-Regel

$$Q^M(f;a,b)=(b-a)f\left(\frac{a+b}{2}\right)$$

Note: Interpoliert nur durch eine Konstante (1 Knoten → exakt für Grad $0$), hat aber trotzdem Ordnung 2, denn sie ist auch für lineare $f$ exakt:

• links von $m$ ist das Rechteck zu hoch, rechts zu tief (oder umgekehrt) → die beiden Dreiecke sind gleich gross und heben sich auf (siehe Bild oben).
• das ist genau Satz 5.3.11 mit $1$ Knoten: Ordnung $\ge 1$ und gerade $\Rightarrow 2$.
• Mittelpunkt = Gauss-Quadratur mit $s=1$ Knoten (Kap. 5.5).

5.4.2 Trapez-Regel

Trapez-Regel

$$Q^T(f;a,b)=\frac{b-a}{2}(f(a)+f(b))$$

Fehler: ${\int }_a^{b}fdx-Q^T(f)=-\frac{1}{12}(b-a{)}^3{f}^{\prime \prime }(\xi )$ für ein $\xi \in [a,b]$.

Wir legen eine Gerade durch $(a,f(a))$ und $(b,f(b))$ und integrieren diese:

$$p(x)=f(a)+\frac{x-a}{b-a}(f(b)-f(a))$$

Note: Das ist “Startwert + Steigung $\cdot$ zurückgelegter Weg”: $\frac{f(b)-f(a)}{b-a}$ ist die Steigung, $x-a$ der Abstand von $a$. Check: $p(a)=f(a)$, $p(b)=f(b)$.

• In der Lagrange-Form: $p(x)=f(a)\underset{{\ell }_0}{\underbrace{\frac{b-x}{b-a}}}+f(b)\underset{{\ell }_1}{\underbrace{\frac{x-a}{b-a}}}$ und ${\int }_a^b{\ell }_0={\int }_a^b{\ell }_1=\frac{b-a}{2}$ → die Gewichte.
• Fläche unter der Geraden = Trapez mit Breite $b-a$ und mittlerer Höhe $\frac{f(a)+f(b)}{2}$.

Note: 2 Knoten, symmetrisch → Ordnung $2$ (Grad $1$ exakt, $x^2$ nicht). Der Fehler ist etwa doppelt so gross wie bei Mittelpunkt und hat das umgekehrte Vorzeichen (Mittelpunkt: $+\frac{1}{24}(b-a{)}^3{f}^{\prime \prime }(\xi )$, Standardresultat, nicht im Skript).

5.4.3 Simpson-Regel

Simpson-Regel

$$Q^S(f;a,b)=\frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right)$$

Fehler: ${\int }_a^{b}fdx-Q^S(f)=-\frac{1}{90}{\left(\frac{b-a}{2}\right)}^5{f}^{(4)}(\xi )$ für ein $\xi \in [a,b]$.

Note: Parabel durch $a,m,b$ → exakt bis Grad $2$, also Ordnung $\ge 3$. Symmetrisch → Ordnung gerade → Ordnung 4, d.h. exakt bis Grad $3$.

• $x^3$ ist ungerade bzgl. $m$ → beide Seiten heben sich auf, die Formel bekommt Grad $3$ gratis.
• $x^4$ wird nicht mehr exakt integriert: auf $[-1,1]$ gibt Simpson $\frac{2}{3}$, das Integral ist $\frac{2}{5}$.

5.4.4 Summierte Formeln

Mit $h=\frac{b-a}{N}$, $x_i=a+ih$:

Summierte Mittelpunkt-, Trapez- und Simpson-Regel

$$M_N(f)=\sum _{i=0}^{N-1}hf\left(\frac{x_i+x_{i+1}}{2}\right)$$ $$T_N(f)=\sum _{i=0}^{N-1}\frac{h}{2}(f(x_i)+f(x_{i+1}))=\frac{h}{2}\left(f(a)+2\sum _{i=1}^{N-1}f(x_i)+f(b)\right)$$ $$S_N(f)=\frac{h}{6}\left(f(a)+2\sum _{i=1}^{N-1}f(x_i)+4\sum _{i=1}^{N}f\left(\frac{x_{i-1}+x_i}{2}\right)+f(b)\right)$$

Note: Innere Knoten gehören zu zwei Teilintervallen → ihr Gewicht verdoppelt sich (Trapez: $\frac{h}{2}\to h$). Simpson braucht $2N+1$ Auswertungen.
Fehler summierte Trapezregel: $E\le \frac{h^2}{12}(b-a){\max }_{x\in [a,b]}|f^{\prime \prime }(x)|$.

Konvergenz der summierten Formeln (Satz 5.4.4)

Basiert $Q$ auf einer lokalen Formel der Ordnung $p$, dann gibt es $C>0$ mit

$$\left|{\int }_{I}f(t)dt-Q(f)\right|\le C{h}^p\Vert f^{(p)}{\Vert }_{L^{\infty }(I)}\text{fu¨r alle }f\in C^p(I)$$

→ Die Ordnung sieht man im log-log-Plot als Steigung. Aber: nur wenn $f$ glatt genug ist! Für $\sqrt{t}$ (Knick bei $0$) kommen alle Formeln nur auf $O(h^{1.5})$.

import numpy as np
 
def trapezoidal(f, a, b, N):
    x, h = np.linspace(a, b, N + 1, retstep=True)
    return h / 2 * (f(x[0]) + 2 * np.sum(f(x[1:-1])) + f(x[-1]))
 
def simpson(f, a, b, N):                       # N Teilintervalle, 2N+1 Punkte
    x, h = np.linspace(a, b, 2 * N + 1, retstep=True)   # h ist hier schon (b-a)/(2N)
    return h / 3 * (f(x[0]) + 4 * np.sum(f(x[1::2])) + 2 * np.sum(f(x[2:-1:2])) + f(x[-1]))

(Skript Code 5.4.6 / 5.4.7)

Note: Überraschung (Bsp. 5.4.8): Für periodische, analytische $f$ über eine ganze Periode konvergiert die summierte Trapezregel exponentiell (nicht nur $O(h^2)$).

• Grund: $T_N$ ist exakt für $e^{2\pi ikt}$ mit $|k|<N$ (das ist dieselbe Rechnung wie bei der DFT, Kap. 3.2).
• nur über eine Teilperiode (z.B. $[0,\frac{1}{2}]$) → wieder algebraisch.

Note: Newton-Cotes höherer Ordnung nicht verwenden (Bem. 5.4.13): Mehr äquidistante Knoten = Interpolation hohen Grades auf äquidistanten Punkten → Runge-Phänomen. Ab 9 Knoten gibt es negative Gewichte → Auslöschung. In der Praxis nur Mittelpunkt, Trapez, Simpson (summiert), sonst Gauss / Clenshaw-Curtis.

5.4.5 Ausblick: Endpunkte als Knoten (Radau, Lobatto)

Mit $s$ frei wählbaren Knoten ist die maximale Ordnung $2s$ (Gauss, Satz 5.5.2). Fixiert man Endpunkte als Knoten, verliert man Freiheitsgrade (Bem. 5.5.11):

	fixe Knoten	max. Ordnung	Beispiel
Gauss	keine	$2s$	Mittelpunkt ($s=1$, Ordnung $2$)
Radau	ein Endpunkt	$2s-1$
Lobatto	$a$ und $b$	$2s-2$	Trapez ($s=2$, Ordnung $2$), Simpson ($s=3$, Ordnung $4$)

→ Trapez und Simpson sind die Lobatto-Formeln mit $2$ bzw. $3$ Knoten, Mittelpunkt ist keine (enthält $a,b$ nicht), sondern Gauss.