In Anwendungen der Mathematik müssen häufig Riemann-Integrale für stückweise stetige Funktionen berechnet werden. In vielen Fällen ist eine geschlossene Lösung eines solchen Integrals nicht bekannt, so dass es näherungsweise numerisch gelöst werden muss. Die numerische Lösung eines Integrals bezeichnet man auch als numerische Quadratur . In diesem Abschnitt sollen eine Reihe von Formeln zur numerischen Integration hergeleitet und untersucht werden.
Das Integral über eine stückweise stetige Funktion kann bekanntlich als Summe von Integralen über stetige Funktionen beschrieben werden, so dass wir uns auf die Betrachtung stetiger Funktionen beschränken können. Dazu definieren wir den Operator
I
:
C
[
a
,
b
]
→
R
{\displaystyle {\mathcal {I}}:C[a,b]\to \mathbb {R} }
mit
I
(
f
)
:=
∫
a
b
f
(
x
)
d
x
{\displaystyle {\mathcal {I}}(f):=\int \limits _{a}^{b}f(x)\,dx}
für
f
∈
C
[
a
,
b
]
{\displaystyle f\in C[a,b]}
. Dieser ist linear, da für alle
f
,
g
∈
C
[
a
,
b
]
{\displaystyle f,g\in C[a,b]}
und
α
,
β
∈
R
{\displaystyle \alpha ,\beta \in \mathbb {R} }
I
(
α
f
+
β
g
)
=
∫
a
b
[
α
f
(
x
)
+
β
g
(
x
)
]
d
x
=
α
∫
a
b
f
(
x
)
d
x
+
β
∫
a
b
g
(
x
)
d
x
=
α
I
(
f
)
+
β
I
(
g
)
{\displaystyle {\mathcal {I}}(\alpha f+\beta g)=\int \limits _{a}^{b}[\alpha f(x)+\beta g(x)]\,dx=\alpha \int \limits _{a}^{b}f(x)\,dx+\beta \int \limits _{a}^{b}g(x)\,dx=\alpha {\mathcal {I}}(f)+\beta {\mathcal {I}}(g)}
gilt und er ist positiv , d. h. man hat
f
∈
C
[
a
,
b
]
,
f
≥
0
⇒
I
(
f
)
≥
0
,
{\displaystyle f\in C[a,b],f\geq 0\Rightarrow {\mathcal {I}}(f)\geq 0,}
wobei
f
≥
0
{\displaystyle f\geq 0}
bedeutet, dass
f
(
x
)
≥
0
,
x
∈
[
a
,
b
]
{\displaystyle f(x)\geq 0,x\in [a,b]}
ist. Gesucht sind nun einfach auszuwertende Formeln, die jedem
f
∈
C
[
a
,
b
]
{\displaystyle f\in C[a,b]}
einen Näherungswert
I
^
(
f
)
{\displaystyle {\hat {\mathcal {I}}}(f)}
für den Wert des Integrals zuordnen und zwar so, dass der Quadraturfehler
I
(
f
)
−
I
^
(
f
)
{\displaystyle {\mathcal {I}}(f)-{\hat {\mathcal {I}}}(f)}
möglichst klein ist.
Unter einer Quadraturformel
I
n
:
C
[
a
,
b
]
→
R
{\displaystyle {\mathcal {I}}_{n}:C[a,b]\to \mathbb {R} }
zur Berechnung des bestimmten Integrals
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
versteht man eine Summe
(8.1)
I
n
(
f
)
:=
(
b
−
a
)
∑
i
=
0
n
σ
i
f
(
x
i
)
{\displaystyle {\mathcal {I}}_{n}(f):=(b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})}
für
f
∈
C
[
a
,
b
]
{\displaystyle f\in C[a,b]}
mit bekannten Gewichten
σ
i
∈
R
{\displaystyle \sigma _{i}\in \mathbb {R} }
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und Stützstellen bzw. Knoten
x
i
∈
[
a
,
b
]
{\displaystyle x_{i}\in [a,b]}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
, wobei
x
i
≠
x
j
{\displaystyle x_{i}\neq x_{j}}
(
i
≠
j
)
{\displaystyle (i\neq j)}
sei.
Wenn wir die Abhängigkeit der Gewichte und Stützstellen von der Wahl von
n
{\displaystyle n}
darstellen wollen, schreiben wir statt
σ
i
{\displaystyle \sigma _{i}}
und
x
i
{\displaystyle x_{i}}
auch
σ
i
(
n
)
{\displaystyle \sigma _{i}^{(n)}}
und
x
i
(
n
)
{\displaystyle x_{i}^{(n)}}
. Wie
I
{\displaystyle {\mathcal {I}}}
ist auch eine Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
ein linearer Operator, denn man hat
I
n
(
α
f
+
β
g
)
=
(
b
−
a
)
∑
i
=
0
n
σ
i
[
α
f
(
x
i
)
+
β
g
(
x
i
)
]
=
α
(
b
−
a
)
∑
i
=
0
n
σ
i
f
(
x
i
)
+
β
(
b
−
a
)
∑
i
=
0
n
σ
i
g
(
x
i
)
{\displaystyle {\mathcal {I}}_{n}(\alpha f+\beta g)=(b-a)\sum _{i=0}^{n}\sigma _{i}[\alpha f(x_{i})+\beta g(x_{i})]=\alpha (b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})+\beta (b-a)\sum _{i=0}^{n}\sigma _{i}g(x_{i})}
(8.2)
=
α
I
n
(
f
)
+
β
I
n
(
g
)
{\displaystyle =\alpha {\mathcal {I}}_{n}(f)+\beta {\mathcal {I}}_{n}(g)}
für alle
f
,
g
∈
C
[
a
,
b
]
{\displaystyle f,g\in C[a,b]}
und
α
,
β
∈
R
{\displaystyle \alpha ,\beta \in \mathbb {R} }
. Sind die Gewichte nichtnegativ, d. h. hat man
σ
i
≥
0
{\displaystyle \sigma _{i}\geq 0}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
, so ist ferner mit
I
{\displaystyle {\mathcal {I}}}
auch
I
n
{\displaystyle {\mathcal {I}}_{n}}
positiv und gilt also
f
∈
C
[
a
,
b
]
,
f
≥
0
⇒
I
n
(
f
)
≥
0.
{\displaystyle f\in C[a,b],f\geq 0\Rightarrow {\mathcal {I}}_{n}(f)\geq 0.}
Wir definieren weiter:
Eine Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
hat mindestens den Genauigkeitsgrad
r
∈
N
0
{\displaystyle r\in \mathbb {N} _{0}}
, wenn
(8.3)
I
n
(
x
j
)
=
I
(
x
j
)
,
j
=
0
,
1
,
…
,
r
{\displaystyle {\mathcal {I}}_{n}(x^{j})={\mathcal {I}}(x^{j}),\quad j=0,1,\ldots ,r}
gilt. Im Fall, dass zusätzlich
I
n
(
x
r
+
1
)
≠
I
(
x
r
+
1
)
{\displaystyle {\mathcal {I}}_{n}(x^{r+1})\neq {\mathcal {I}}(x^{r+1})}
richtig ist, sagt man, dass
I
n
{\displaystyle {\mathcal {I}}_{n}}
den Genauigkeitsgrad
r
∈
N
0
{\displaystyle r\in \mathbb {N} _{0}}
hat.
Da
I
{\displaystyle {\mathcal {I}}}
und
I
n
{\displaystyle {\mathcal {I}}_{n}}
lineare Operatoren sind, folgt aus (8.3)
(8.4)
I
n
(
∑
j
=
0
r
a
j
x
j
)
=
∑
j
=
0
r
a
j
I
n
(
x
j
)
=
∑
j
=
0
r
a
j
I
(
x
j
)
=
I
(
∑
j
=
0
r
a
j
x
j
)
{\displaystyle {\mathcal {I}}_{n}\left(\sum _{j=0}^{r}a_{j}x^{j}\right)=\sum _{j=0}^{r}a_{j}{\mathcal {I}}_{n}(x^{j})=\sum _{j=0}^{r}a_{j}{\mathcal {I}}(x^{j})={\mathcal {I}}\left(\sum _{j=0}^{r}a_{j}x^{j}\right)}
für alle
a
j
∈
R
{\displaystyle a_{j}\in \mathbb {R} }
(
j
=
0
,
1
,
…
,
r
)
{\displaystyle (j=0,1,\ldots ,r)}
, und damit der folgende Satz, wobei
Π
r
{\displaystyle \Pi _{r}}
wieder den Raum aller Polynome vom Grad
≤
r
{\displaystyle \leq r}
bezeichnet:
I
n
{\displaystyle {\mathcal {I}}_{n}}
ist genau dann eine Quadraturformel von mindestens dem Genauigkeitsgrad
r
{\displaystyle r}
, wenn gilt:
I
n
(
p
)
=
I
(
p
)
,
p
∈
Π
r
.
{\displaystyle {\mathcal {I}}_{n}(p)={\mathcal {I}}(p),\quad p\in \Pi _{r}.}
Wir bemerken in diesem Zusammenhang:
Zu
n
+
1
{\displaystyle n+1}
Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
mit
x
i
≠
x
k
{\displaystyle x_{i}\neq x_{k}}
(
i
≠
k
)
{\displaystyle (i\neq k)}
gibt es genau eine Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
, welche mindestens den Genauigkeitsgrad
n
∈
N
0
{\displaystyle n\in \mathbb {N} _{0}}
hat, d. h. für die gilt:
(8.5)
I
n
(
p
)
=
I
(
p
)
,
p
∈
Π
n
.
{\displaystyle {\mathcal {I}}_{n}(p)={\mathcal {I}}(p),\quad p\in \Pi _{n}.}
Diese hat die Gewichte
(8.6)
σ
i
:=
1
b
−
a
I
(
L
i
)
,
i
=
0
,
1
,
…
,
n
,
{\displaystyle \sigma _{i}:={\frac {1}{b-a}}{\mathcal {I}}(L_{i}),\quad i=0,1,\ldots ,n,}
wobei
L
i
∈
Π
n
{\displaystyle L_{i}\in \Pi _{n}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
die zu den
x
i
{\displaystyle x_{i}}
gehörenden Lagrangeschen Basispolynome sind (vgl. Definition 6.2).
Für die durch die Stützstellen
x
i
{\displaystyle x_{i}}
und Gewichte
σ
i
{\displaystyle \sigma _{i}}
in (8.6) definierte Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
gilt für
k
=
0
,
1
,
…
,
n
{\displaystyle k=0,1,\ldots ,n}
(8.7)
I
n
(
L
k
)
=
(
b
−
a
)
∑
i
=
0
n
σ
i
L
k
(
x
i
)
=
(
b
−
a
)
∑
i
=
0
n
σ
i
δ
k
i
=
(
b
−
a
)
σ
k
=
I
(
L
k
)
.
{\displaystyle {\mathcal {I}}_{n}(L_{k})=(b-a)\sum _{i=0}^{n}\sigma _{i}L_{k}(x_{i})=(b-a)\sum _{i=0}^{n}\sigma _{i}\delta _{ki}=(b-a)\sigma _{k}={\mathcal {I}}(L_{k}).}
Da sich jedes Polynom vom Grad
≤
n
{\displaystyle \leq n}
auf eindeutige Weise als Linearkombination der
L
k
{\displaystyle L_{k}}
(
k
=
0
,
1
,
…
,
n
)
{\displaystyle (k=0,1,\ldots ,n)}
darstellen lässt (vgl. (6.6)) und
I
{\displaystyle {\mathcal {I}}}
sowie
I
n
{\displaystyle {\mathcal {I}}_{n}}
lineare Operatoren sind, folgt damit analog zu (8.4) die Beziehung (8.5). Die so definierte Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
ist eindeutig. Denn für jede andere Quadraturformel
J
n
{\displaystyle {\mathcal {J}}_{n}}
mit Gewichten
σ
^
k
{\displaystyle {\hat {\sigma }}_{k}}
und Genauigkeitsgrad
n
≥
0
{\displaystyle n\geq 0}
hat man wegen
L
k
∈
Π
n
{\displaystyle L_{k}\in \Pi _{n}}
die Identität
J
n
(
L
k
)
=
I
(
L
k
)
{\displaystyle {\mathcal {J}}_{n}(L_{k})={\mathcal {I}}(L_{k})}
und bekommt man analog zu (8.7)
J
n
(
L
k
)
=
(
b
−
a
)
σ
^
k
{\displaystyle {\mathcal {J}}_{n}(L_{k})=(b-a){\hat {\sigma }}_{k}}
, so dass
(
b
−
a
)
σ
^
k
=
I
(
L
k
)
{\displaystyle (b-a){\hat {\sigma }}_{k}={\mathcal {I}}(L_{k})}
und demnach
σ
^
k
=
σ
k
{\displaystyle {\hat {\sigma }}_{k}=\sigma _{k}}
folgt.
q.e.d.
Weiter stellen wir fest:
Ist
I
n
{\displaystyle {\mathcal {I}}_{n}}
eine Quadraturformel, die einen Genauigkeitsgrad
r
∈
N
0
{\displaystyle r\in \mathbb {N} _{0}}
hat, so folgt für ihre Gewichte
σ
i
{\displaystyle \sigma _{i}}
∑
i
=
0
n
σ
i
=
1.
{\displaystyle \sum _{i=0}^{n}\sigma _{i}=1.}
Da
I
n
{\displaystyle {\mathcal {I}}_{n}}
einen Genauigkeitsgrad
r
≥
0
{\displaystyle r\geq 0}
hat, folgt
(
b
−
a
)
∑
i
=
0
n
σ
i
=
I
n
(
1
)
=
I
(
1
)
=
∫
a
b
1
d
x
=
b
−
a
.
{\displaystyle (b-a)\sum _{i=0}^{n}\sigma _{i}={\mathcal {I}}_{n}(1)={\mathcal {I}}(1)=\int \limits _{a}^{b}1\,dx=b-a.}
q.e.d.
Bezüglich der Konvergenz der durch eine Quadraturformel
I
n
(
f
)
:=
(
b
−
a
)
∑
i
=
0
n
σ
i
(
n
)
f
(
x
i
(
n
)
{\displaystyle {\mathcal {I}}_{n}(f):=(b-a)\sum _{i=0}^{n}\sigma _{i}^{(n)}f(x_{i}^{(n)}}
erzeugten Näherungswerte
I
n
(
f
)
{\displaystyle {\mathcal {I}}_{n}(f)}
gegen den exakten Wert des Integrals
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
für
n
→
∞
{\displaystyle n\to \infty }
kann man allgemein den folgenden Satz angeben, den wir hier jedoch nicht beweisen können (für einen Beweis siehe H. Heuser: Funktionalanalysis, Teubner, Stuttgart, 1992, S. 268).
Man hat
lim
n
→
∞
I
n
(
f
)
=
I
(
f
)
,
f
∈
C
[
a
,
b
]
{\displaystyle \lim _{n\to \infty }{\mathcal {I}}_{n}(f)={\mathcal {I}}(f),\quad f\in C[a,b]}
genau dann, wenn gilt:
(a)
∑
i
=
1
n
|
σ
i
(
n
)
|
≤
M
,
n
∈
N
,
{\displaystyle \sum _{i=1}^{n}\left|\sigma _{i}^{(n)}\right|\leq M,\quad n\in \mathbb {N} ,}
(b)
lim
n
→
∞
I
n
(
x
j
)
=
I
(
x
j
)
,
j
∈
N
0
.
{\displaystyle \lim _{n\to \infty }{\mathcal {I}}_{n}(x^{j})={\mathcal {I}}(x^{j}),\quad j\in \mathbb {N} _{0}.}
Mit Hilfe von Satz 8.5 erschließt man ferner:
Es sei
I
n
{\displaystyle {\mathcal {I}}_{n}}
eine Quadraturformel mit Gewichten
σ
i
(
n
)
≥
0
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle \sigma _{i}^{(n)}\geq 0\ (i=0,1,\ldots ,n)}
für alle
n
∈
N
{\displaystyle n\in \mathbb {N} }
und einem Genauigkeitsgrad
n
≥
0
{\displaystyle n\geq 0}
. Dann hat man
lim
n
→
∞
I
n
(
f
)
=
I
(
f
)
,
f
∈
C
[
a
,
b
]
{\displaystyle \lim _{n\to \infty }{\mathcal {I}}_{n}(f)={\mathcal {I}}(f),\quad f\in C[a,b]}
genau dann, wenn gilt:
lim
n
→
∞
I
n
(
x
j
)
=
I
(
x
j
)
,
j
∈
N
0
.
{\displaystyle \lim _{n\to \infty }{\mathcal {I}}_{n}(x^{j})={\mathcal {I}}(x^{j}),\quad j\in \mathbb {N} _{0}.}
Es seien nun
f
∈
C
[
a
,
b
]
{\displaystyle f\in C[a,b]}
und
x
i
∈
[
a
,
b
]
,
i
=
0
,
1
,
…
,
n
{\displaystyle x_{i}\in [a,b],i=0,1,\ldots ,n}
mit
x
i
≠
x
k
{\displaystyle x_{i}\neq x_{k}}
(
i
≠
k
)
{\displaystyle (i\neq k)}
gegeben und
Q
n
∈
Π
n
{\displaystyle Q_{n}\in \Pi _{n}}
bezeichne das (eindeutige) Interpolationspolynom zu den Stützpunkten
(
x
i
,
f
(
x
i
)
)
,
i
=
0
,
1
,
…
,
n
{\displaystyle (x_{i},f(x_{i})),i=0,1,\ldots ,n}
. Sind
L
i
∈
Π
n
{\displaystyle L_{i}\in \Pi _{n}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
wieder die zu den
n
+
1
{\displaystyle n+1}
Stützstellen
x
i
{\displaystyle x_{i}}
gehörenden Lagrangeschen Basispolynome, so kann das Interpolationspolynom
Q
n
{\displaystyle Q_{n}}
damit gemäß (6.7)in der Form
Q
n
(
x
)
=
∑
i
=
0
n
f
(
x
i
)
L
i
(
x
)
{\displaystyle Q_{n}(x)=\sum _{i=0}^{n}f(x_{i})L_{i}(x)}
geschrieben werden. Wir definieren nun:
Eine Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
mit
I
n
(
f
)
:=
I
(
Q
n
)
=
∫
a
b
Q
n
(
x
)
d
x
=
(
b
−
a
)
∑
i
=
0
n
{
1
b
−
a
∫
a
b
L
i
(
x
)
d
x
}
f
(
x
i
)
,
{\displaystyle {\mathcal {I}}_{n}(f):={\mathcal {I}}(Q_{n})=\int \limits _{a}^{b}Q_{n}(x)\,dx=(b-a)\sum _{i=0}^{n}\left\{{\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(x)\,dx\right\}f(x_{i}),}
d. h. mit Gewichten
(8.8)
σ
i
:=
1
b
−
a
I
(
L
i
)
,
i
=
0
,
1
,
…
,
n
{\displaystyle \sigma _{i}:={\frac {1}{b-a}}{\mathcal {I}}(L_{i}),\quad i=0,1,\ldots ,n}
heißt interpolatorische Quadraturformel .
Wegen der Übereinstimmung der Gewichte in (8.8) und (8.6) können wir mit Satz 8.4 schließen:
Eine interpolatorische Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
hat mindestens den Genauigkeitsgrad
n
∈
N
0
{\displaystyle n\in \mathbb {N} _{0}}
und ist zu den gegebenen Stützstellen die einzige Quadraturformel mit einem Genauigkeitsgrad
≥
n
{\displaystyle \geq n}
.
Ferner können wir zeigen:
Eine interpolatorische Quadraturformel In besitzt die Gestalt
I
n
(
f
)
=
(
b
−
a
)
∑
i
=
0
n
σ
i
f
(
x
i
)
{\displaystyle {\mathcal {I}}_{n}(f)=(b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})}
mit
(8.9)
σ
i
:=
∫
0
1
∏
k
=
0
k
≠
i
n
t
−
t
k
t
i
−
t
k
d
t
{\displaystyle \sigma _{i}:=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {t-t_{k}}{t_{i}-t_{k}}}\,dt}
mit
t
k
:=
x
k
−
a
b
−
a
.
{\displaystyle t_{k}:={\frac {x_{k}-a}{b-a}}.}
Mit
t
k
{\displaystyle t_{k}}
wie in (8.9) lassen sich die Gewichte
σ
i
{\displaystyle \sigma _{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
aus (8.8) mit Hilfe der Substitution
x
:=
(
b
−
a
)
t
+
a
{\displaystyle x:=(b-a)t+a}
umschreiben in
σ
i
=
1
b
−
a
I
(
L
i
)
=
1
b
−
a
∫
a
b
L
i
(
x
)
d
x
=
1
b
−
a
∫
a
b
∏
k
=
0
k
≠
i
n
x
−
x
k
x
i
−
x
k
d
x
{\displaystyle \sigma _{i}={\frac {1}{b-a}}{\mathcal {I}}(L_{i})={\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}\prod _{k=0 \atop k\neq i}^{n}{\frac {x-x_{k}}{x_{i}-x_{k}}}\,dx}
(8.10)
=
∫
0
1
∏
k
=
0
k
≠
i
n
(
b
−
a
)
t
−
(
b
−
a
)
t
k
(
b
−
a
)
t
i
−
(
b
−
a
)
t
k
d
t
=
∫
0
1
∏
k
=
0
k
≠
i
n
t
−
t
k
t
i
−
t
k
d
t
.
{\displaystyle =\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {(b-a)t-(b-a)t_{k}}{(b-a)t_{i}-(b-a)t_{k}}}\,dt=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {t-t_{k}}{t_{i}-t_{k}}}\,dt.}
q.e.d.
Die Transformation von
x
{\displaystyle x}
nach
t
{\displaystyle t}
in (8.9) ist sinnvoll, da damit die Gewichte
σ
i
{\displaystyle \sigma _{i}}
in einer interpolatorischen Quadraturformel von den Intervallgrenzen
a
{\displaystyle a}
und
b
{\displaystyle b}
unabhängig werden und nur von der relativen Verteilung der Stützstellen in
[
a
,
b
]
{\displaystyle [a,b]}
abhängen.
Wir wollen nun auf spezielle interpolatorische Quadraturformeln, die Newton-Cotes-Formeln , eingehen. Diese ergeben sich durch äquidistante Wahl der Stützstellen in
[
a
,
b
]
{\displaystyle [a,b]}
. Insbesondere erhält man die abgeschlossenen Newton-Cotes-Formeln , wenn die Randpunkte des Intervalls
[
a
,
b
]
{\displaystyle [a,b]}
selbst Stützstellen sind, wenn also für
n
≥
1
g
i
l
t
::
(
8.11
)
<
m
a
t
h
>
x
i
:=
a
+
i
h
(
i
=
0.1
,
…
,
n
)
,
h
:=
(
b
−
a
)
/
n
.
{\displaystyle n\geq 1gilt::(8.11)<math>x_{i}:=a+ih\quad (i=0.1,\ldots ,n),\qquad h:=(b-a)/n.}
Bei den offenen Newton-Cotes-Formeln sind die Randpunkte von
[
a
,
b
]
{\displaystyle [a,b]}
selbst keine Stützstellen, so dass man
x
i
:=
a
+
(
i
+
1
)
h
(
i
=
0
,
1
,
…
,
n
)
,
h
:=
(
b
−
a
)
/
(
n
+
2
)
{\displaystyle x_{i}:=a+(i+1)h\quad (i=0,1,\ldots ,n),\qquad h:=(b-a)/(n+2)}
hat. Wir wollen hier nur die abgeschlossenen Newton-Cotes-Formeln genauer untersuchen.
Für die Gewichte
σ
i
{\displaystyle \sigma _{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
der abgeschlossenen Newton-Cotes-Formeln gilt:
(8.12)
σ
n
−
i
=
σ
i
=
1
n
∫
0
n
∏
k
=
0
k
≠
i
n
s
−
k
i
−
k
d
s
.
{\displaystyle \sigma _{n-i}=\sigma _{i}={\frac {1}{n}}\int \limits _{0}^{n}\prod _{k=0 \atop k\neq i}^{n}{\frac {s-k}{i-k}}\,ds.}
Die zweite Identität in (8.12) folgt mit
t
k
:=
x
k
−
a
b
−
a
=
k
h
b
−
a
=
k
n
{\displaystyle t_{k}:={\frac {x_{k}-a}{b-a}}={\frac {kh}{b-a}}={\frac {k}{n}}}
aus (8.9) mit der Substitution
t
:=
s
/
n
{\displaystyle t:=s/n}
, denn man hat
σ
i
=
∫
0
1
∏
k
=
0
k
≠
i
n
t
−
k
n
i
n
−
k
n
d
t
=
1
n
∫
0
n
∏
k
=
0
k
≠
i
n
s
−
k
i
−
k
d
s
.
{\displaystyle \sigma _{i}=\int \limits _{0}^{1}\prod _{k=0 \atop k\neq i}^{n}{\frac {t-{\frac {k}{n}}}{{\frac {i}{n}}-{\frac {k}{n}}}}\,dt={\frac {1}{n}}\int \limits _{0}^{n}\prod _{k=0 \atop k\neq i}^{n}{\frac {s-k}{i-k}}\,ds.}
Somit müssen wir noch die erste Identität in (8.12) zeigen. Dazu sei
i
∈
{
0
,
…
,
n
}
{\displaystyle i\in \{0,\ldots ,n\}}
. Sind
L
i
∈
Π
n
{\displaystyle L_{i}\in \Pi _{n}}
die Lagrangeschen Basispolynome, so ist
L
n
−
i
∈
Π
n
{\displaystyle L_{n-i}\in \Pi _{n}}
und
Q
(
x
)
:=
L
i
(
b
+
a
−
x
)
∈
Π
n
{\displaystyle Q(x):=L_{i}(b+a-x)\in \Pi _{n}}
sowie
Q
(
x
n
−
j
)
=
L
i
(
b
+
a
−
x
n
−
j
)
=
L
i
(
b
+
a
−
[
a
+
(
n
−
j
)
b
−
a
n
]
)
=
L
i
(
a
+
j
b
−
a
n
)
=
L
i
(
x
j
)
{\displaystyle Q(x_{n-j})=L_{i}(b+a-x_{n-j})=L_{i}\left(b+a-\left[a+(n-j){\frac {b-a}{n}}\right]\right)=L_{i}\left(a+j{\frac {b-a}{n}}\right)=L_{i}(x_{j})}
=
δ
i
j
=
δ
n
−
i
,
n
−
j
=
L
n
−
i
(
x
n
−
j
)
{\displaystyle =\delta _{ij}=\delta _{n-i,n-j}=L_{n-i}(x_{n-j})}
für
j
=
0
,
1
,
…
,
n
{\displaystyle j=0,1,\ldots ,n}
. Da
L
n
−
i
{\displaystyle L_{n-i}}
und
Q
{\displaystyle Q}
demnach offenbar Interpolationspolynome zu den Punkten
(
x
n
−
j
,
δ
n
−
i
,
n
−
j
)
,
j
=
0
,
1
,
…
,
n
{\displaystyle (x_{n-j},\delta _{n-i,n-j}),j=0,1,\ldots ,n}
sind, muss wegen der Eindeutigkeit des Interpolationspolynoms
L
n
−
i
≡
Q
{\displaystyle L_{n-i}\equiv Q}
gelten, so dass wir schließlich mit der Substitution
t
:=
b
+
a
−
x
{\displaystyle t:=b+a-x}
Folgendes erhalten (vgl. (8.10)):
σ
n
−
i
=
1
b
−
a
∫
a
b
L
n
−
i
(
x
)
d
x
=
1
b
−
a
∫
a
b
Q
(
x
)
d
x
=
1
b
−
a
∫
a
b
L
i
(
b
+
a
−
x
)
d
x
=
1
b
−
a
∫
a
b
L
i
(
t
)
d
t
=
σ
i
.
{\displaystyle \sigma _{n-i}={\frac {1}{b-a}}\int \limits _{a}^{b}L_{n-i}(x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}Q(x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(b+a-x)\,dx={\frac {1}{b-a}}\int \limits _{a}^{b}L_{i}(t)\,dt=\sigma _{i}.}
q.e.d.
Wir geben nun einige Spezialfälle der abgeschlossenen Newton-Cotes-Formeln an.
(1) Für
n
=
1
{\displaystyle n=1}
hat eine interpolatorische Quadraturformel nach Satz 8.10 die Gestalt
I
1
(
f
)
=
(
b
−
a
)
[
σ
0
f
(
x
0
)
+
σ
1
f
(
x
1
)
]
.
{\displaystyle {\mathcal {I}}_{1}(f)=(b-a)[\sigma _{0}f(x_{0})+\sigma _{1}f(x_{1})].}
Dabei ergeben sich für die zugehörige abgeschlossene Newton-Cotes-Formel mit
h
=
b
−
a
{\displaystyle h=b-a}
die Stützstellen
x
0
=
a
{\displaystyle x_{0}=a}
und
x
1
=
b
{\displaystyle x_{1}=b}
und wegen
σ
1
=
σ
0
{\displaystyle \sigma _{1}=\sigma _{0}}
(Lemma 8.11) und
σ
0
+
σ
1
=
1
{\displaystyle \sigma _{0}+\sigma _{1}=1}
(Satz 8.5) die Gewichte
σ
1
=
σ
0
=
1
2
.
{\displaystyle \sigma _{1}=\sigma _{0}={\frac {1}{2}}.}
Man erhält so die (Sehnen-) Trapezregel
(8.13)
I
1
(
f
)
:=
b
−
a
2
[
f
(
a
)
+
f
(
b
)
]
.
{\displaystyle {\mathcal {I}}_{1}(f):={\frac {b-a}{2}}[f(a)+f(b)].}
(2) Für
n
=
2
{\displaystyle n=2}
hat man mit
h
=
(
b
−
a
)
/
2
{\displaystyle h=(b-a)/2}
die Stützstellen
x
0
=
a
,
x
1
=
(
a
+
b
)
/
2
{\displaystyle x_{0}=a,x_{1}=(a+b)/2}
und
x
2
=
b
{\displaystyle x_{2}=b}
und unter Verwendung von Lemma 8.11 und anschließend Satz 8.5 die Gewichte
σ
2
=
σ
0
=
1
2
0
2
(
s
−
1
)
(
s
−
2
)
(
0
−
1
)
(
0
−
2
)
d
s
=
1
4
0
2
(
s
2
−
3
s
+
2
)
d
s
=
1
4
[
8
3
−
6
+
4
]
=
1
6
,
{\displaystyle \sigma _{2}=\sigma _{0}={\frac {1}{2}}_{0}^{2}{\frac {(s-1)(s-2)}{(0-1)(0-2)}}\,ds={\frac {1}{4}}_{0}^{2}(s^{2}-3s+2)\,ds={\frac {1}{4}}\left[{\frac {8}{3}}-6+4\right]={\frac {1}{6}},}
σ
1
=
1
−
σ
0
−
σ
2
=
1
−
2
6
=
2
3
.
{\displaystyle \sigma _{1}=1-\sigma _{0}-\sigma _{2}=1-{\frac {2}{6}}={\frac {2}{3}}.}
Unter Verwendung von Satz 8.10 ergibt sich so die Simpson-Regel bzw. Keplersche Fassregel
(8.14)
I
2
(
f
)
:=
b
−
a
6
[
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
]
{\displaystyle {\mathcal {I}}_{2}(f):={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]}
(3) Der Fall
n
=
3
{\displaystyle n=3}
führt auf die Newtonsche 3/8-Regel
I
3
(
f
)
:=
b
−
a
8
[
f
(
a
)
+
3
f
(
2
a
+
b
3
)
+
3
f
(
a
+
2
b
3
)
+
f
(
b
)
]
.
{\displaystyle {\mathcal {I}}_{3}(f):={\frac {b-a}{8}}\left[f(a)+3f\left({\frac {2a+b}{3}}\right)+3f\left({\frac {a+2b}{3}}\right)+f(b)\right].}
(4) Für
n
=
4
{\displaystyle n=4}
bekommt man die Milne-Regel
I
4
(
f
)
:=
b
−
a
90
[
7
f
(
a
)
+
32
f
(
3
a
+
b
4
)
+
12
f
(
2
a
+
2
b
4
)
+
32
f
(
a
+
3
b
4
)
+
7
f
(
b
)
]
.
{\displaystyle {\mathcal {I}}_{4}(f):={\frac {b-a}{90}}\left[7f(a)+32f\left({\frac {3a+b}{4}}\right)+12f\left({\frac {2a+2b}{4}}\right)+32f\left({\frac {a+3b}{4}}\right)+7f(b)\right].}
Als Beispiel berechnen wir ein Integral näherungsweise mit der Simpson-Regel.
Es seien
f
(
x
)
:=
1
/
(
1
+
x
2
)
,
a
=
0
{\displaystyle f(x):=1/(1+x^{2}),a=0}
und
b
=
1
{\displaystyle b=1}
, so dass
I
(
f
)
=
∫
0
1
1
1
+
x
2
d
x
{\displaystyle {\mathcal {I}}(f)=\int \limits _{0}^{1}{\frac {1}{1+x^{2}}}\,dx}
ist. Die Simpson-Regel liefert dafür den Näherungswert
I
2
(
f
)
=
1
6
[
f
(
0
)
+
4
f
(
1
2
)
+
f
(
1
)
]
=
1
6
(
1
+
4
⋅
4
5
+
1
2
)
=
47
60
=
0.783
33.
{\displaystyle {\mathcal {I}}_{2}(f)={\frac {1}{6}}\left[f(0)+4f\left({\frac {1}{2}}\right)+f(1)\right]={\frac {1}{6}}\left(1+4\cdot {\frac {4}{5}}+{\frac {1}{2}}\right)={\frac {47}{60}}=0.783\,33.}
Der exakte Wert des Integrals lautet hier
I
(
f
)
=
arctan
(
x
)
|
x
=
0
x
=
1
=
arctan
(
1
)
=
0.785
40.
{\displaystyle {\mathcal {I}}(f)=\arctan(x)|_{x=0}^{x=1}=\arctan(1)=0.785\,40.}
Für
n
≤
7
{\displaystyle n\leq 7}
sind die Gewichte in den abgeschlossenen Newton-Cotes-Formeln nichtnegativ und sind diese Quadraturformeln demzufolge positiv. Für
n
=
8
{\displaystyle n=8}
und
n
≥
10
{\displaystyle n\geq 10}
treten negative Gewichte auf und ist damit als Folge von Satz 8.5
∑
i
=
0
n
|
σ
i
(
n
)
|
>
1
,
{\displaystyle \sum _{i=0}^{n}\left|\sigma _{i}^{(n)}\right|>1,}
was zu einer Verstärkung von Rundungsfehlern bei den Funktionswerten
f
(
x
i
)
{\displaystyle f(x_{i})}
führt. Die Verwendung der abgeschlossenen Newton-Cotes-Formeln für
n
≥
8
{\displaystyle n\geq 8}
ist daher nicht zu empfehlen. Für die (abgeschlossenen) Newton-Cotes-Formeln lässt sich sogar
lim
n
→
∞
∑
i
=
0
n
|
σ
i
(
n
)
|
=
1
{\displaystyle \lim _{n\to \infty }\sum _{i=0}^{n}\left|\sigma _{i}^{(n)}\right|=1}
beweisen (Satz von Kusmin), so dass man aus dem Satz 8.6 von Szegö die Existenz eines
f
∈
C
[
a
,
b
]
{\displaystyle f\in C[a,b]}
schließen kann, für das die Konvergenz
lim
n
→
∞
I
n
(
f
)
=
I
(
f
)
{\displaystyle \lim _{n\to \infty }{\mathcal {I}}_{n}(f)={\mathcal {I}}(f)}
nicht gilt. Letzteres lässt ja auch der Satz 6.24 von Faber generell für interpolatorische Quadraturformeln vermuten. Eine Erhöhung von
n
{\displaystyle n}
bei den (abgeschlossenen) Newton-Cotes-Formeln muss also nicht zwangsläufig zu einer genaueren Näherung
I
n
(
f
)
{\displaystyle {\mathcal {I}}_{n}(f)}
von
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
führen.
Wir geben hier noch einige weitere interpolatorische Quadraturformeln an.
(1) Für
n
=
0
{\displaystyle n=0}
und
x
0
:=
a
{\displaystyle x_{0}:=a}
oder
x
0
:=
b
{\displaystyle x_{0}:=b}
muss wegen Satz 8.5
σ
0
=
1
{\displaystyle \sigma _{0}=1}
gelten, so dass man alternativ folgende beiden Rechteckregeln erhält:
(8.15)
I
0
(
f
)
:=
(
b
−
a
)
f
(
a
)
,
I
0
(
f
)
:=
(
b
−
a
)
f
(
b
)
.
{\displaystyle {\mathcal {I}}_{0}(f):=(b-a)f(a),\quad {\mathcal {I}}_{0}(f):=(b-a)f(b).}
(2) Für
n
=
0
{\displaystyle n=0}
bekommt man im Fall der offenen Newton-Cotes-Formeln
h
:=
(
b
−
a
)
/
2
{\displaystyle h:=(b-a)/2}
und
x
0
:=
a
+
b
−
a
2
=
a
+
b
2
,
σ
0
=
∑
i
=
0
n
σ
i
=
1
{\displaystyle x_{0}:=a+{\frac {b-a}{2}}={\frac {a+b}{2}},\quad \sigma _{0}=\sum _{i=0}^{n}\sigma _{i}=1}
und damit eine weitere Rechteckregel, die Mittelpunktregel
I
0
(
f
)
:=
(
b
−
a
)
f
(
a
+
b
2
)
.
{\displaystyle {\mathcal {I}}_{0}(f):=(b-a)f\left({\frac {a+b}{2}}\right).}
(3) Die offene Newton-Cotes-Formel für
n
=
1
{\displaystyle n=1}
lautet mit
h
:=
(
b
−
a
)
/
3
{\displaystyle h:=(b-a)/3}
,
x
0
:=
a
+
1
3
(
b
−
a
)
,
x
1
:=
a
+
2
3
(
b
−
a
)
{\displaystyle x_{0}:=a+{\frac {1}{3}}(b-a),\quad x_{1}:=a+{\frac {2}{3}}(b-a)}
und den Gewichten
σ
0
=
σ
1
=
1
/
2
{\displaystyle \sigma _{0}=\sigma _{1}=1/2}
, die man aus der Formel (8.9) errechnet, wie folgt:
I
1
(
f
)
:=
b
−
a
2
[
f
(
2
a
+
b
3
)
+
f
(
a
+
2
b
3
)
]
.
{\displaystyle {\mathcal {I}}_{1}(f):={\frac {b-a}{2}}\left[f\left({\frac {2a+b}{3}}\right)+f\left({\frac {a+2b}{3}}\right)\right].}
(4) Die offene Newton-Cotes-Formel für
n
=
2
{\displaystyle n=2}
lautet mit
h
:=
(
b
−
a
)
/
4
{\displaystyle h:=(b-a)/4}
,
x
0
:=
a
+
1
4
(
b
−
a
)
,
x
1
:=
a
+
1
2
(
b
−
a
)
,
x
2
:=
a
+
3
4
(
b
−
a
)
{\displaystyle x_{0}:=a+{\frac {1}{4}}(b-a),\quad x_{1}:=a+{\frac {1}{2}}(b-a),\quad x_{2}:=a+{\frac {3}{4}}(b-a)}
und den mit Hilfe von (8.9) zu berechnenden Gewichten wie folgt:
I
2
(
f
)
:=
b
−
a
3
[
2
f
(
3
a
+
b
4
)
−
f
(
a
+
b
2
)
+
2
f
(
3
a
+
b
4
)
]
.
{\displaystyle {\mathcal {I}}_{2}(f):={\frac {b-a}{3}}\left[2f\left({\frac {3a+b}{4}}\right)-f\left({\frac {a+b}{2}}\right)+2f\left({\frac {3a+b}{4}}\right)\right].}
Man beachte, dass sie ein negatives Gewicht beinhaltet.
8.2.3 Quadraturfehler und Genauigkeitsgrad
Bearbeiten
Für den durch eine beliebige interpolatorische Quadraturformel in Bezug auf den exakten Wert des Integrals entstehenden Fehler, kann man die im folgenden Satz angegebene Abschätzung beweisen.
Es sei
I
n
{\displaystyle {\mathcal {I}}_{n}}
eine interpolatorische Quadraturformel mit Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
, welche mindestens den Genauigkeitsgrad
r
≥
n
{\displaystyle r\geq n}
besitze und es sei
f
∈
C
(
r
+
1
)
[
a
,
b
]
{\displaystyle f\in C^{(r+1)}[a,b]}
. Dann gilt
(8.16)
|
I
(
f
)
−
I
n
(
f
)
|
≤
γ
r
(
b
−
a
)
r
+
2
(
r
+
1
)
!
max
ξ
∈
[
a
,
b
]
|
f
(
r
+
1
)
(
ξ
)
|
{\displaystyle |{\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)|\leq \gamma _{r}{\frac {(b-a)^{r+2}}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|}
für
γ
r
:=
min
t
n
+
1
,
…
,
t
r
∈
[
0
,
1
]
∫
0
1
∏
k
=
0
r
|
t
−
t
k
|
d
t
{\displaystyle \gamma _{r}:=\min _{t_{n+1},\ldots ,t_{r}\in [0,1]}\int \limits _{0}^{1}\prod _{k=0}^{r}|t-t_{k}|\,dt}
mit
(8.17)
t
k
:=
x
k
−
a
b
−
a
(
k
=
0
,
1
,
…
,
n
)
.
{\displaystyle t_{k}:={\frac {x_{k}-a}{b-a}}\quad (k=0,1,\ldots ,n).}
Hat man insbesondere für die
t
k
{\displaystyle t_{k}}
(
k
=
0
,
…
,
n
)
{\displaystyle (k=0,\ldots ,n)}
aus (8.17) und frei wählbare
t
k
{\displaystyle t_{k}}
(
k
=
n
+
1
,
…
,
r
)
{\displaystyle (k=n+1,\ldots ,r)}
mit
s
(
t
)
:=
∏
k
=
0
r
(
t
−
t
k
)
{\displaystyle s(t):=\prod _{k=0}^{r}(t-t_{k})}
die Beziehung
s
(
t
)
≥
0
,
t
∈
[
0
,
1
]
{\displaystyle s(t)\geq 0,t\in [0,1]}
oder
s
(
t
)
≤
0
,
t
∈
[
0
,
1
]
{\displaystyle s(t)\leq 0,t\in [0,1]}
, dann folgt mit
γ
^
r
:=
∫
0
1
s
(
t
)
d
t
{\displaystyle {\hat {\gamma }}_{r}:=\int \limits _{0}^{1}s(t)\,dt}
und einem
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
die Fehlerdarstellung
(8.18)
I
(
f
)
−
I
n
(
f
)
=
γ
^
r
(
b
−
a
)
r
+
2
(
r
+
1
)
!
f
(
r
+
1
)
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\hat {\gamma }}_{r}{\frac {(b-a)^{r+2}}{(r+1)!}}f^{(r+1)}(\xi ).}
Seien
x
i
∈
[
a
,
b
]
,
i
=
n
+
1
,
…
,
r
{\displaystyle x_{i}\in [a,b],i=n+1,\ldots ,r}
zunächst beliebig gewählt, so dass die
x
i
{\displaystyle x_{i}}
(
i
=
0
,
…
,
n
,
n
+
1
,
…
,
r
)
{\displaystyle (i=0,\ldots ,n,n+1,\ldots ,r)}
paarweise verschieden sind und sei
Q
r
∈
Π
r
{\displaystyle Q_{r}\in \Pi _{r}}
das Interpolationspolynom zur den Stützpunkten
(
x
i
,
f
(
x
i
)
)
,
i
=
0
,
1
,
…
,
r
{\displaystyle (x_{i},f(x_{i})),i=0,1,\ldots ,r}
. Da
I
n
{\displaystyle {\mathcal {I}}_{n}}
den Genauigkeitsgrad
r
{\displaystyle r}
hat, gilt dann
I
n
(
f
)
=
(
b
−
a
)
∑
i
=
0
n
σ
i
f
(
x
i
)
=
(
b
−
a
)
∑
i
=
0
n
σ
i
Q
r
(
x
i
)
=
I
n
(
Q
r
)
=
I
(
Q
r
)
{\displaystyle {\mathcal {I}}_{n}(f)=(b-a)\sum _{i=0}^{n}\sigma _{i}f(x_{i})=(b-a)\sum _{i=0}^{n}\sigma _{i}Q_{r}(x_{i})={\mathcal {I}}_{n}(Q_{r})=I(Qr)}
und demnach
I
(
f
)
−
I
n
(
f
)
=
I
(
f
)
−
I
(
Q
r
)
=
∫
a
b
[
f
(
x
)
−
Q
r
(
x
)
]
d
x
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\mathcal {I}}(f)-I(Q_{r})=\int \limits _{a}^{b}[f(x)-Q_{r}(x)]\,dx.}
Mit
ω
(
x
)
:=
(
x
−
x
0
)
⋯
(
x
−
x
n
)
,
φ
(
x
)
:=
(
x
−
x
n
+
1
)
⋯
(
x
−
x
r
)
{\displaystyle \omega (x):=(x-x_{0})\cdots (x-x_{n}),\quad \varphi (x):=(x-x_{n+1})\cdots (x-x_{r})}
und unter Verwendung von Satz 6.11 hat man für ein
ξ
(
x
)
∈
[
a
,
b
]
{\displaystyle \xi (x)\in [a,b]}
f
(
x
)
−
Q
r
(
x
)
=
1
(
r
+
1
)
!
ω
(
x
)
φ
(
x
)
f
(
r
+
1
)
(
ξ
(
x
)
)
.
{\displaystyle f(x)-Q_{r}(x)={\frac {1}{(r+1)!}}\omega (x)\varphi (x)f^{(r+1)}(\xi (x)).}
Da die linke Seite der letzten Gleichung stetig in
x
{\displaystyle x}
ist, ist es auch die rechte Seite und darum hat man
(8.19)
I
(
f
)
−
I
n
(
f
)
=
1
(
r
+
1
)
!
∫
a
b
[
ω
(
x
)
φ
(
x
)
f
(
r
+
1
)
(
ξ
(
x
)
)
]
d
x
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\frac {1}{(r+1)!}}\int \limits _{a}^{b}\left[\omega (x)\varphi (x)f^{(r+1)}(\xi (x))\right]\,dx.}
Man beachte nun, dass die
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
durch die Quadraturregel festgelegt sind. Wir wollen abschließend zeigen, dass für die Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
n
+
1
,
…
,
r
)
{\displaystyle (i=n+1,\ldots ,r)}
die anfangs gemachte Voraussetzung hinsichtlich der paarweisen Unterschiedlichkeit fallen gelassen werden kann. Es seien daher jetzt letztere Punkte vollkommen beliebig aus [a, b] gewählt. Für jedes
m
∈
N
{\displaystyle m\in \mathbb {N} }
können wir dann Punkte
x
i
(
m
)
{\displaystyle x_{i}^{(m)}}
(
i
=
n
+
1
,
…
,
r
)
{\displaystyle (i=n+1,\ldots ,r)}
finden, die zusammen mit den
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
paarweise unterschiedlich sind und für die
lim
m
→
∞
x
i
(
m
)
=
x
i
(
i
=
n
+
1
,
…
,
r
)
{\displaystyle \lim _{m\to \infty }x_{i}^{(m)}=x_{i}\quad (i=n+1,\ldots ,r)}
gilt. Setzen wir
φ
m
(
x
)
:=
(
x
−
x
n
+
1
(
m
)
)
⋯
(
x
−
x
r
(
m
)
)
,
{\displaystyle \varphi _{m}(x):=(x-x_{n+1}^{(m)})\cdots (x-x_{r}^{(m)}),}
so hat man unter Verwendung des ersten Teils des Beweises
|
I
(
f
)
−
I
n
(
f
)
|
≤
1
(
r
+
1
)
!
max
ξ
∈
[
a
,
b
]
|
f
(
r
+
1
)
(
ξ
)
|
∫
a
b
|
ω
(
x
)
φ
m
(
x
)
|
d
x
{\displaystyle |{\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)|\leq {\frac {1}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|\int \limits _{a}^{b}|\omega (x)\varphi _{m}(x)|\,dx}
≤
1
(
r
+
1
)
!
max
ξ
∈
[
a
,
b
]
|
f
(
r
+
1
)
(
ξ
)
|
{
∫
a
b
|
ω
(
x
)
φ
(
x
)
|
d
x
+
∫
a
b
|
ω
(
x
)
|
|
φ
m
(
x
)
−
φ
(
x
)
|
d
x
⏟
→
0
(
m
→
∞
)
}
{\displaystyle \leq {\frac {1}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|\{\int \limits _{a}^{b}|\omega (x)\varphi (x)|\,dx+\underbrace {\int \limits _{a}^{b}|\omega (x)||\varphi _{m}(x)-\varphi (x)|\,dx} _{\to 0\quad (m\to \infty )}\}}
(8.20)
≤
1
(
r
+
1
)
!
max
ξ
∈
[
a
,
b
]
|
f
(
r
+
1
)
(
ξ
)
|
∫
a
b
|
ω
(
x
)
φ
(
x
)
|
d
x
.
{\displaystyle \leq {\frac {1}{(r+1)!}}\max _{\xi \in [a,b]}\left|f^{(r+1)}(\xi )\right|\int \limits _{a}^{b}|\omega (x)\varphi (x)|\,dx.}
Wählt man nun
x
i
{\displaystyle x_{i}}
(
i
=
n
+
1
,
…
,
r
)
{\displaystyle (i=n+1,\ldots ,r)}
so, dass der Wert des letzten Integrals minimal wird und wendet man die Substitution
x
:=
(
b
−
a
)
t
+
a
{\displaystyle x:=(b-a)t+a}
an, so gelangt man schließlich zu
γ
r
:=
min
x
n
+
1
,
…
,
x
r
∈
[
a
,
b
]
∫
a
b
∏
i
=
0
r
|
x
−
x
i
|
d
x
=
(
b
−
a
)
r
+
2
min
t
n
+
1
,
…
,
t
r
∈
[
0
,
1
]
∫
0
1
∏
i
=
0
r
|
t
−
t
i
|
d
t
,
{\displaystyle \gamma _{r}:=\min _{x_{n+1},\ldots ,x_{r}\in [a,b]}\int \limits _{a}^{b}\prod _{i=0}^{r}|x-x_{i}|\,dx=(b-a)^{r+2}\min _{t_{n+1},\ldots ,t_{r}\in [0,1]}\int \limits _{0}^{1}\prod _{i=0}^{r}|t-t_{i}|\,dt,}
womit die Abschätzung (8.16) gezeigt ist.
Ist nun mit gewissen Punkten
x
i
{\displaystyle x_{i}}
(
i
=
n
+
1
,
…
,
r
)
{\displaystyle (i=n+1,\ldots ,r)}
(8.21)
ω
(
x
)
φ
(
x
)
≥
0
,
x
∈
[
a
,
b
]
,
{\displaystyle \omega (x)\varphi (x)\geq 0,\quad x\in [a,b],}
so erhält man aus (8.20)
I
(
f
)
−
I
n
(
f
)
≤
1
(
r
+
1
)
!
max
η
∈
[
a
,
b
]
|
f
(
r
+
1
)
(
η
)
|
∫
a
b
[
ω
(
x
)
φ
(
x
)
]
d
x
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)\leq {\frac {1}{(r+1)!}}\max _{\eta \in [a,b]}\left|f^{(r+1)}(\eta )\right|\int \limits _{a}^{b}[\omega (x)\varphi (x)]\,dx.}
Weiter gewinnt man mit (8.19)
I
(
f
)
−
I
n
(
f
)
≥
1
(
r
+
1
)
!
min
η
∈
[
a
,
b
]
|
f
(
r
+
1
)
(
η
)
|
∫
a
b
[
ω
(
x
)
φ
(
x
)
]
d
x
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)\geq {\frac {1}{(r+1)!}}\min _{\eta \in [a,b]}\left|f^{(r+1)}(\eta )\right|\int \limits _{a}^{b}[\omega (x)\varphi (x)]\,dx.}
Der Zwischenwertsatz, angewandt auf die Funktion
f
(
r
+
1
)
{\displaystyle f^{(r+1)}}
, liefert somit für ein
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
I
(
f
)
−
I
n
(
f
)
=
1
(
r
+
1
)
!
f
(
r
+
1
)
(
ξ
)
∫
a
b
[
ω
(
x
)
φ
(
x
)
]
d
x
,
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\frac {1}{(r+1)!}}f^{(r+1)}(\xi )\int \limits _{a}^{b}[\omega (x)\varphi (x)]\,dx,}
so dass die Substitution
x
:=
(
b
−
a
)
t
+
a
{\displaystyle x:=(b-a)t+a}
in diesem Fall zu der Formel (8.18) führt. Analog schließt man im Fall, dass „
≤
{\displaystyle \leq }
“ statt „
≥
{\displaystyle \geq }
“ in (8.21) vorliegt.
q.e.d.
Wir nutzen im Folgenden aus, dass nach Korollar 8.9 der Genauigkeitsgrad einer interpolatorischen Quadraturformel
I
n
m
i
n
d
e
s
t
e
n
s
<
m
a
t
h
>
r
:=
n
i
s
t
.
(
1
)
S
e
i
<
m
a
t
h
>
f
∈
C
1
[
a
,
b
]
{\displaystyle {\mathcal {I}}_{n}mindestens<math>r:=nist.(1)Sei<math>f\in C^{1}[a,b]}
und
I
0
(
f
)
:=
(
b
−
a
)
f
(
a
)
{\displaystyle {\mathcal {I}}_{0}(f):=(b-a)f(a)}
die Rechteckregel aus (8.15). Aus (8.18) gewinnt man für
I
0
{\displaystyle {\mathcal {I}}_{0}}
mit
r
=
n
=
0
{\displaystyle r=n=0}
, mit
x
0
:=
a
{\displaystyle x_{0}:=a}
bzw.
t
0
:=
0
{\displaystyle t_{0}:=0}
sowie mit
∏
k
=
0
0
(
t
−
t
k
)
=
t
≥
0
,
t
∈
[
0
,
1
]
{\displaystyle \prod _{k=0}^{0}(t-t_{k})=t\geq 0,\quad t\in [0,1]}
und
γ
^
0
:=
∫
0
1
t
d
t
=
1
2
{\displaystyle {\hat {\gamma }}_{0}:=\int \limits _{0}^{1}t\,dt={\frac {1}{2}}}
die Fehlerdarstellung
I
(
f
)
−
I
0
(
f
)
=
(
b
−
a
)
2
2
f
′
(
ξ
)
,
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{0}(f)={\frac {(b-a)^{2}}{2}}f'(\xi ),}
wobei
ξ
{\displaystyle \xi }
ein Punkt aus
[
a
,
b
]
{\displaystyle [a,b]}
ist. Entsprechend erhält man für die Rechteckregel
I
0
(
f
)
:=
(
b
−
a
)
f
(
b
)
{\displaystyle {\mathcal {I}}_{0}(f):=(b-a)f(b)}
mit
r
=
n
=
0
,
x
0
:=
b
{\displaystyle r=n=0,x_{0}:=b}
bzw.
t
0
:=
1
{\displaystyle t_{0}:=1}
sowie mit
∏
k
=
0
0
(
t
−
t
k
)
=
t
−
1
≤
0
,
t
∈
[
0
,
1
]
{\displaystyle \prod _{k=0}^{0}(t-t_{k})=t-1\leq 0,\quad t\in [0,1]}
und
γ
^
0
:=
∫
0
1
(
t
−
1
)
d
t
=
−
1
2
{\displaystyle {\hat {\gamma }}_{0}:=\int \limits _{0}^{1}(t-1)\,dt=-{\frac {1}{2}}}
die Fehlerdarstellung
I
(
f
)
−
I
0
(
f
)
=
−
(
b
−
a
)
2
2
f
′
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{0}(f)=-{\frac {(b-a)^{2}}{2}}f'(\xi ).}
(2) Im Fall der Trapezregel
I
1
(
f
)
:=
b
−
a
2
[
f
(
a
)
+
f
(
b
)
]
{\displaystyle {\mathcal {I}}_{1}(f):={\frac {b-a}{2}}[f(a)+f(b)]}
gilt für
f
∈
C
2
[
a
,
b
]
{\displaystyle f\in C^{2}[a,b]}
mit einem
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
die Fehlerdarstellung
I
(
f
)
−
I
1
(
f
)
=
−
(
b
−
a
)
3
12
f
″
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{1}(f)=-{\frac {(b-a)^{3}}{12}}f''(\xi ).}
Denn mit
r
=
n
=
1
,
x
0
:=
a
,
x
1
:=
b
{\displaystyle r=n=1,x_{0}:=a,x_{1}:=b}
bzw.
t
0
:=
0
,
t
1
:=
1
{\displaystyle t_{0}:=0,t_{1}:=1}
hat man
∏
k
=
0
1
(
t
−
t
k
)
=
t
(
t
−
1
)
≤
0
,
t
∈
[
0
,
1
]
{\displaystyle \prod _{k=0}^{1}(t-t_{k})=t(t-1)\leq 0,\quad t\in [0,1]}
sowie
γ
^
1
:=
∫
0
1
t
(
t
−
1
)
d
t
=
1
3
−
1
2
=
−
1
6
.
{\displaystyle {\hat {\gamma }}_{1}:=\int \limits _{0}^{1}t(t-1)\,dt={\frac {1}{3}}-{\frac {1}{2}}=-{\frac {1}{6}}.}
Der Genauigkeitsgrad einer interpolatorischen Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
ist mindestens
r
:=
n
{\displaystyle r:=n}
. Für gerade
n
{\displaystyle n}
hat man im Fall der abgeschlossenen Newton-Cotes-Formeln sogar das folgende Resultat (für den Beweis siehe Plato, S. 103):
Die abgeschlossene Newton-Cotes-Formel
I
n
{\displaystyle {\mathcal {I}}_{n}}
besitzt für gerades
n
≥
2
{\displaystyle n\geq 2}
den (exakten) Genauigkeitsgrad
r
:=
n
+
1
{\displaystyle r:=n+1}
.
Letzteres Ergebnis können wir z. B. für die Fehlerdarstellung der Simpson-Regel verwenden.
Es sei
f
∈
C
4
[
a
,
b
]
{\displaystyle f\in C^{4}[a,b]}
. Dann hat man für
n
=
2
{\displaystyle n=2}
und
r
=
3
{\displaystyle r=3}
mit
x
0
:=
a
,
x
1
:=
(
a
+
b
)
/
2
,
x
2
:=
b
{\displaystyle x_{0}:=a,x_{1}:=(a+b)/2,x_{2}:=b}
bzw.
t
0
:=
0
,
t
1
:=
1
/
2
,
t
2
:=
1
{\displaystyle t_{0}:=0,t_{1}:=1/2,t_{2}:=1}
und mit dem gewählten Punkt
t
3
:=
1
/
2
{\displaystyle t_{3}:=1/2}
∏
k
=
0
3
(
t
−
t
k
)
=
t
(
t
−
1
2
)
2
(
t
−
1
)
≤
0
,
t
∈
[
0
,
1
]
{\displaystyle \prod _{k=0}^{3}(t-t_{k})=t\left(t-{\frac {1}{2}}\right)^{2}(t-1)\leq 0,\quad t\in [0,1]}
sowie
γ
^
3
:=
∫
0
1
[
t
(
t
−
1
2
)
2
(
t
−
1
)
]
d
t
=
−
1
120
.
{\displaystyle {\hat {\gamma }}_{3}:=\int \limits _{0}^{1}\left[t\left(t-{\frac {1}{2}}\right)^{2}(t-1)\right]\,dt=-{\frac {1}{120}}.}
Also ergibt sich für die Simpson-Regel
I
2
(
f
)
:=
b
−
a
6
[
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
]
{\displaystyle {\mathcal {I}}_{2}(f):={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]}
mit einem
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
der Quadraturfehler
I
(
f
)
−
I
2
(
f
)
=
−
(
b
−
a
)
5
4
!
⋅
120
f
(
4
)
(
ξ
)
=
−
(
b
−
a
)
5
2880
f
(
4
)
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{2}(f)=-{\frac {(b-a)^{5}}{4!\cdot 120}}f^{(4)}(\xi )=-{\frac {(b-a)^{5}}{2880}}f^{(4)}(\xi ).}
Wie bereits in Abschnitt 8.2.2 erläutert wurde, garantiert eine Erhöhung von
n
{\displaystyle n}
keineswegs, dass die Newton-Cotes-Formeln Näherungswerte zunehmender Genauigkeit für
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
liefern. Um Letzteres zu erreichen, müssen wir daher anders vorgehen. Und zwar teilen wir zunächst das Intervall
[
a
,
b
]
{\displaystyle [a,b]}
mittels Stützstellen
x
k
:=
a
+
k
h
(
k
=
0
,
1
,
…
,
N
)
,
h
:=
b
−
a
N
{\displaystyle x_{k}:=a+kh\quad (k=0,1,\ldots ,N),\qquad h:={\frac {b-a}{N}}}
in
N
{\displaystyle N}
gleiche Stücke auf, so dass sich insbesondere
h
=
x
k
+
1
−
x
k
{\displaystyle h=x_{k+1}-x_{k}}
für alle
k
∈
{
0
,
1
,
…
,
N
−
1
}
{\displaystyle k\in \{0,1,\ldots ,N-1\}}
ergibt. Dann nähern wir das Integral
I
(
f
)
=
∫
a
b
f
(
x
)
d
x
=
∑
k
=
0
N
−
1
∫
x
k
x
k
+
1
f
(
x
)
d
x
{\displaystyle {\mathcal {I}}(f)=\int \limits _{a}^{b}f(x)\,dx=\sum _{k=0}^{N-1}\int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx}
durch
J
n
(
f
)
:=
∑
k
=
0
N
−
1
∫
x
k
x
k
+
1
Q
n
(
x
)
d
x
{\displaystyle {\mathcal {J}}_{n}(f):=\sum _{k=0}^{N-1}\int \limits _{x_{k}}^{x_{k+1}}Q_{n}(x)\,dx}
an, wobei
Q
n
∈
Π
n
{\displaystyle Q_{n}\in \Pi _{n}}
das (eindeutige) Interpolationspolynom zu
n
+
1
{\displaystyle n+1}
paarweise verschiedenen, in jedem Intervall
[
x
k
,
x
k
+
1
]
{\displaystyle [x_{k},x_{k+1}]}
in gleichen Abständen gewählten Stützpunkten ist (vgl. Definition 8.8). Wir wählen also eine interpolatorische Quadraturformel und ersetzen jedes der Integrale
∫
x
k
x
k
+
1
f
(
x
)
d
x
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx}
durch den sich damit ergebenden Wert. Eine so gewonnene Quadraturformel bezeichnet man als summierte Quadraturformel . Wir wollen solche Formeln nun genauer betrachten, wobei wir uns hier auf die abgeschlossenen Newton-Cotes-Formeln zu deren Generierung beschränken wollen. Letztere Wahl legt die Stützpunkte in jedem Intervall
[
x
k
,
x
k
+
1
]
{\displaystyle [x_{k},x_{k+1}]}
durch (8.11) fest, wobei dort
a
:=
x
k
{\displaystyle a:=x_{k}}
und
b
:=
x
k
+
1
{\displaystyle b:=x_{k+1}}
zu wählen ist.
Wir beginnen mit den beiden Rechteckregeln aus (8.15). Für diese erhält man
∫
x
k
x
k
+
1
Q
0
(
x
)
d
x
=
h
f
(
x
k
)
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}Q_{0}(x)\,dx=hf(x_{k})}
bzw.
∫
x
k
x
k
+
1
Q
0
(
x
)
d
x
=
h
f
(
x
k
+
1
)
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}Q_{0}(x)\,dx=hf(x_{k+1})}
so dass Summation über
k
{\displaystyle k}
die folgenden summierten Rechteckregeln liefert:
J
0
(
h
)
:=
h
∑
k
=
0
N
−
1
f
(
x
k
)
,
J
^
0
(
h
)
:=
h
∑
k
=
0
N
−
1
f
(
x
k
+
1
)
.
{\displaystyle {\mathcal {J}}_{0}(h):=h\sum _{k=0}^{N-1}f(x_{k}),\quad {\hat {\mathcal {J}}}_{0}(h):=h\sum _{k=0}^{N-1}f(x_{k+1}).}
Für diese gelten die nachstehenden Fehlerabschätzungen.
Es sei
f
∈
C
1
[
a
,
b
]
{\displaystyle f\in C^{1}[a,b]}
. Dann gibt es
ξ
,
ξ
^
∈
[
a
,
b
]
{\displaystyle \xi ,{\hat {\xi }}\in [a,b]}
, so dass gilt:
(8.22)
I
(
f
)
−
J
0
(
h
)
=
b
−
a
2
h
f
′
(
ξ
)
,
I
(
f
)
−
J
^
0
(
h
)
=
b
−
a
2
h
f
′
(
ξ
^
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {J}}_{0}(h)={\frac {b-a}{2}}hf'(\xi ),\quad {\mathcal {I}}(f)-{\hat {\mathcal {J}}}_{0}(h)={\frac {b-a}{2}}hf'({\hat {\xi }}).}
Aus Beispiel 8.16 (1) ergibt sich für
k
=
0
,
1
,
…
,
N
−
1
{\displaystyle k=0,1,\ldots ,N-1}
die Existenz eines
ξ
k
∈
[
a
,
b
]
{\displaystyle \xi _{k}\in [a,b]}
mit
∫
x
k
x
k
+
1
f
(
x
)
d
x
−
h
f
(
x
k
)
=
h
2
2
f
′
(
ξ
k
)
.
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx-hf(x_{k})={\frac {h^{2}}{2}}f'(\xi _{k}).}
Summation über
k
{\displaystyle k}
führt auf
I
(
f
)
−
J
0
(
h
)
=
h
2
2
∑
k
=
0
N
−
1
f
′
(
ξ
k
)
=
b
−
a
2
h
1
N
∑
k
=
0
N
−
1
f
′
(
ξ
k
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {J}}_{0}(h)={\frac {h^{2}}{2}}\sum _{k=0}^{N-1}f'(\xi _{k})={\frac {b-a}{2}}h{\frac {1}{N}}\sum _{k=0}^{N-1}f'(\xi _{k}).}
Aufgrund von
N
min
x
∈
[
a
,
b
]
f
′
(
x
)
≤
∑
k
=
0
N
−
1
f
′
(
ξ
k
)
≤
N
max
x
∈
[
a
,
b
]
f
′
(
x
)
{\displaystyle N\min _{x\in [a,b]}f'(x)\leq \sum _{k=0}^{N-1}f'(\xi _{k})\leq N\max _{x\in [a,b]}f'(x)}
bzw.
min
x
∈
[
a
,
b
]
f
′
(
x
)
≤
1
N
∑
k
=
0
N
−
1
f
′
(
ξ
k
)
≤
max
x
∈
[
a
,
b
]
f
′
(
x
)
{\displaystyle \min _{x\in [a,b]}f'(x)\leq {\frac {1}{N}}\sum _{k=0}^{N-1}f'(\xi _{k})\leq \max _{x\in [a,b]}f'(x)}
existiert nach dem Zwischenwertsatz ein
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
mit
f
′
(
ξ
)
=
1
N
∑
k
=
0
N
−
1
f
′
(
ξ
k
)
,
{\displaystyle f'(\xi )={\frac {1}{N}}\sum _{k=0}^{N-1}f'(\xi _{k}),}
so dass die erste Fehlerdarstellung in (8.22) folgt. Die zweite zeigt man analog.
q.e.d.
Im Fall der Trapezregel (8.13) hat man
∫
x
k
x
k
+
1
Q
1
(
x
)
d
x
=
h
2
[
f
(
x
k
)
+
f
(
x
k
+
1
)
]
.
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}Q_{1}(x)\,dx={\frac {h}{2}}[f(x_{k})+f(x_{k+1})].}
Summation über
k
{\displaystyle k}
führt auf die summierte Trapezregel
J
1
(
h
)
:=
h
2
(
f
(
a
)
+
2
∑
k
=
1
N
−
1
f
(
x
k
)
+
f
(
b
)
)
{\displaystyle {\mathcal {J}}_{1}(h):={\frac {h}{2}}\left(f(a)+2\sum _{k=1}^{N-1}f(x_{k})+f(b)\right)}
mit der im folgenden Satz angegebenen Fehlerdarstellung.
Es sei
f
∈
C
2
[
a
,
b
]
{\displaystyle f\in C^{2}[a,b]}
. Dann existiert ein
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
mit
I
(
f
)
−
J
1
(
h
)
=
−
b
−
a
12
h
2
f
″
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {J}}_{1}(h)=-{\frac {b-a}{12}}h^{2}f''(\xi ).}
Der Beweis verläuft analog zu dem von Satz 8.19. Nach Beispiel 8.16 (2) gibt es für
k
=
0
,
1
,
…
,
N
−
1
{\displaystyle k=0,1,\ldots ,N-1}
ein
ξ
k
∈
[
a
,
b
]
{\displaystyle \xi _{k}\in [a,b]}
mit
∫
x
k
x
k
+
1
f
(
x
)
d
x
−
h
2
[
f
(
x
k
)
+
f
(
x
k
+
1
)
]
=
−
h
3
12
f
″
(
ξ
k
)
.
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx-{\frac {h}{2}}[f(x_{k})+f(x_{k+1})]=-{\frac {h^{3}}{12}}f''(\xi _{k}).}
Summation über
k
{\displaystyle k}
liefert mit einem
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
I
(
f
)
−
J
1
(
h
)
=
−
b
−
a
12
h
2
1
N
∑
k
=
0
N
−
1
f
″
(
ξ
k
)
=
−
b
−
a
12
h
2
f
″
(
ξ
)
,
{\displaystyle {\mathcal {I}}(f)-{\mathcal {J}}_{1}(h)=-{\frac {b-a}{12}}h^{2}{\frac {1}{N}}\sum _{k=0}^{N-1}f''(\xi _{k})=-{\frac {b-a}{12}}h^{2}f''(\xi ),}
wobei die Existenz eines solchen
ξ
{\displaystyle \xi }
aus der Anwendung des Zwischenwertsatzes auf
f
″
{\displaystyle f''}
geschlossen werden kann.
q.e.d.
Schließlich betrachten wir noch die summierte Simpson-Regel, wobei wir die Darstellung
x
k
:=
a
+
k
h
{\displaystyle x_{k}:=a+kh}
mit
h
:=
(
b
−
a
)
/
N
{\displaystyle h:=(b-a)/N}
für jedes
k
≥
0
{\displaystyle k\geq 0}
verwenden, so dass insbesondere
x
k
+
x
k
+
1
2
=
x
k
+
1
2
(
x
k
+
1
−
x
k
)
=
a
+
k
h
+
1
2
h
=
x
k
+
1
/
2
,
k
=
0
,
1
,
…
,
N
−
1
{\displaystyle {\frac {x_{k}+x_{k+1}}{2}}=x_{k}+{\frac {1}{2}}(x_{k+1}-x_{k})=a+kh+{\frac {1}{2}}h=x_{k+1/2},\quad k=0,1,\ldots ,N-1}
folgt. Die Simpson-Regel, angewandt auf das Intervall
[
x
k
,
x
k
+
1
]
{\displaystyle [x_{k},x_{k+1}]}
, lässt sich somit in der Form
∫
x
k
x
k
+
1
Q
2
(
x
)
=
h
6
[
f
(
x
k
)
+
4
f
(
x
k
+
1
/
2
)
+
f
(
x
k
+
1
)
]
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}Q_{2}(x)={\frac {h}{6}}\left[f(x_{k})+4f(x_{k+1/2})+f(x_{k+1})\right]}
schreiben. Summation über
k
{\displaystyle k}
führt auf die summierte Simpson-Regel
J
2
(
h
)
:=
h
6
(
f
(
a
)
+
4
∑
k
=
0
N
−
1
f
(
x
k
+
1
/
2
)
+
2
∑
k
=
1
N
−
1
f
(
x
k
)
+
f
(
b
)
)
{\displaystyle {\mathcal {J}}_{2}(h):={\frac {h}{6}}\left(f(a)+4\sum _{k=0}^{N-1}f(x_{k+1/2})+2\sum _{k=1}^{N-1}f(x_{k})+f(b)\right)}
Für diese hat man die im folgenden Satz angegebene Fehlerdarstellung.
Es sei
f
∈
C
4
[
a
,
b
]
{\displaystyle f\in C^{4}[a,b]}
. Dann existiert ein
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
mit
I
(
f
)
−
J
2
(
h
)
=
−
b
−
a
2880
h
4
f
(
4
)
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {J}}_{2}(h)=-{\frac {b-a}{2880}}h^{4}f^{(4)}(\xi ).}
Der Beweis verläuft wiederum analog zu dem von Satz 8.19. Nach Beispiel 8.18 gibt es für
k
=
0
,
1
,
…
,
N
−
1
{\displaystyle k=0,1,\ldots ,N-1}
ein
ξ
k
∈
[
a
,
b
]
{\displaystyle \xi _{k}\in [a,b]}
mit
∫
x
k
x
k
+
1
f
(
x
)
d
x
−
h
6
[
f
(
x
k
)
+
4
f
(
x
k
+
1
/
2
)
+
f
(
x
k
+
1
)
]
=
−
h
5
2880
f
(
4
)
(
ξ
k
)
.
{\displaystyle \int \limits _{x_{k}}^{x_{k+1}}f(x)\,dx-{\frac {h}{6}}\left[f(x_{k})+4f(x_{k+1/2})+f(x_{k+1})\right]=-{\frac {h^{5}}{2880}}f^{(4)}(\xi _{k}).}
Summation über
k
{\displaystyle k}
liefert mit einem
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
I
(
f
)
−
J
2
(
h
)
=
−
b
−
a
2880
h
4
1
N
∑
k
=
0
N
−
1
f
(
4
)
(
ξ
k
)
=
−
b
−
a
2880
h
4
f
(
4
)
(
ξ
)
,
{\displaystyle {\mathcal {I}}(f)-{\mathcal {J}}_{2}(h)=-{\frac {b-a}{2880}}h^{4}{\frac {1}{N}}\sum _{k=0}^{N-1}f^{(4)}(\xi _{k})=-{\frac {b-a}{2880}}h^{4}f^{(4)}(\xi ),}
wobei die Existenz von
ξ
{\displaystyle \xi }
aus der Anwendung des Zwischenwertsatzes auf
f
(
4
)
{\displaystyle f^{(4)}}
folgt.
q.e.d.
Zur Auswertung der summierten Rechteckregeln müssen
N
{\displaystyle N}
, für die der summierten Trapezregel
N
+
1
{\displaystyle N+1}
und für die der summierte Simpson-Regel
2
N
+
1
{\displaystyle 2N+1}
Funktionswerte bestimmt werden. Der Rechenaufwand bei Verwendung der summierten Simpson-Regel ist damit etwa doppelt so hoch wie der bei Verwendung einer der drei anderen Regeln. Dennoch ist die summierte Simpson-Regel diesen für hinreichend glatte Funktionen wegen der höheren Fehlerordnung in h vorzuziehen. Denn der Quadraturfehler verhält sich bei ihr wie
O
(
h
4
)
{\displaystyle {\mathcal {O}}(h^{4})}
, während er bei den summierten Rechteckregeln und der summierten Trapezregel proportional zu
h
{\displaystyle h}
bzw.
h
2
{\displaystyle h^{2}}
abnimmt.
Da man die in der jeweiligen Fehlerformel vorkommende Ableitung durch das Maximum des Betrages dieser Ableitung bezüglich aller
x
∈
[
a
,
b
]
{\displaystyle x\in [a,b]}
nach oben abschätzen kann, implizieren die angegebenen Fehlerdarstellungen insbesondere, dass die hier angegebenen summierten Quadraturformeln für
h
→
0
{\displaystyle h\to 0}
gegen den exakten Wert des Integrals
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
konvergieren, wobei mit „
h
→
0
{\displaystyle h\to 0}
“ hier „
h
k
=
(
b
−
a
)
/
N
k
{\displaystyle h_{k}=(b-a)/N_{k}}
mit
N
k
∈
N
{\displaystyle N_{k}\in \mathbb {N} }
und
N
k
→
∞
{\displaystyle N_{k}\to \infty }
“ gemeint ist.
Wir greifen abschließend nochmals Beispiel 8.13 auf.
Es seien wieder
f
(
x
)
:=
1
/
(
1
+
x
2
)
,
a
=
0
{\displaystyle f(x):=1/(1+x^{2}),a=0}
und
b
=
1
{\displaystyle b=1}
, so dass ein Näherungswert für das Integral
I
(
f
)
=
∫
0
1
1
1
+
x
2
d
x
{\displaystyle {\mathcal {I}}(f)=\int \limits _{0}^{1}{\frac {1}{1+x^{2}}}\,dx}
gesucht ist. Weiter wählen wir
N
=
3
{\displaystyle N=3}
und somit
h
=
1
/
3
{\displaystyle h=1/3}
. Der exakte Wert des Integrals lautet
arctan
(
1
)
=
0.785
398
1
{\displaystyle \arctan(1)=0.785\,398\,1}
. Mit der summierten Simpson-Regel ergibt sich der Wert
J
2
(
h
)
=
1
18
[
f
(
0
)
+
4
(
f
(
1
/
6
)
+
f
(
3
/
6
)
+
f
(
5
/
6
)
)
+
2
(
f
(
1
/
3
)
+
f
(
2
/
3
)
)
+
f
(
1
)
]
{\displaystyle {\mathcal {J}}_{2}(h)={\frac {1}{18}}\left[f(0)+4\left(f(1/6)+f(3/6)+f(5/6)\right)+2\left(f(1/3)+f(2/3)\right)+f(1)\right]}
=
1
18
[
1
+
4
(
36
37
+
36
45
+
36
61
)
+
2
(
9
10
+
9
13
)
+
1
2
]
=
0.785
397
94.
{\displaystyle ={\frac {1}{18}}\left[1+4\left({\frac {36}{37}}+{\frac {36}{45}}+{\frac {36}{61}}\right)+2\left({\frac {9}{10}}+{\frac {9}{13}}\right)+{\frac {1}{2}}\right]=0.785\,397\,94.}
Für die summierte Trapezregel
J
1
(
h
)
{\displaystyle {\mathcal {J}}_{1}(h)}
gibt der folgende Satz eine asymptotische Entwicklung nach Potenzen von
h
2
{\displaystyle h^{2}}
an, welche dazu genutzt werden soll, aus einer endlichen Zahl von Auswertungen der summierten Trapezregel eine im Hinblick auf diese Werte genauere Näherung des Integrals
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
zu berechnen. (Der Satz ist z. B. bei Plato bewiesen.)
Für ein
r
≥
0
{\displaystyle r\geq 0}
sei
f
∈
C
2
r
+
2
[
a
,
b
]
{\displaystyle f\in C^{2r+2}[a,b]}
. Die summierte Trapezregel
J
1
(
h
)
:=
h
2
(
f
(
a
)
+
2
∑
k
=
1
N
−
1
f
(
x
k
)
+
f
(
b
)
)
{\displaystyle {\mathcal {J}}_{1}(h):={\frac {h}{2}}\left(f(a)+2\sum _{k=1}^{N-1}f(x_{k})+f(b)\right)}
mit
h
:=
(
b
−
a
)
/
N
{\displaystyle h:=(b-a)/N}
für ein
N
∈
N
{\displaystyle N\in \mathbb {N} }
besitzt die asymptotische Entwicklung
(8.23)
J
1
(
h
)
=
α
0
+
α
1
h
2
+
α
2
h
2
⋅
2
+
…
+
α
r
h
2
r
+
O
(
h
2
r
+
2
)
{\displaystyle {\mathcal {J}}_{1}(h)=\alpha _{0}+\alpha _{1}h^{2}+\alpha _{2}h^{2\cdot 2}+\ldots +\alpha _{r}h^{2r}+{\mathcal {O}}(h^{2r+2})}
für
h
→
0
{\displaystyle h\to 0}
mit
α
0
:=
I
(
f
)
{\displaystyle \alpha _{0}:=I(f)}
und gewissen Koeffizienten
α
i
∈
R
{\displaystyle \alpha _{i}\in \mathbb {R} }
(
i
=
1
,
…
,
r
)
{\displaystyle (i=1,\ldots ,r)}
.
Für periodische Funktionen mit Periode
b
−
a
{\displaystyle b-a}
kann man sogar zeigen, dass
α
i
=
0
{\displaystyle \alpha _{i}=0}
(
i
=
1
,
…
,
r
)
{\displaystyle (i=1,\ldots ,r)}
gilt. In einem solchen Fall kann mit dem in diesem Abschnitt beschriebenen Verfahren keine Verbesserung erzielt werden.
Man beachte, dass man
J
1
(
h
)
{\displaystyle {\mathcal {J}}_{1}(h)}
nur für
h
>
0
{\displaystyle h>0}
mit
h
:=
(
b
−
a
)
/
N
{\displaystyle h:=(b-a)/N}
für eine natürliche Zahl
N
{\displaystyle N}
auswerten kann. Aufgrund von (8.23) (wie auch wegen Satz 8.20) gilt ferner
lim
h
→
0
J
1
(
h
)
=
I
(
f
)
,
{\displaystyle \lim _{h\to 0}{\mathcal {J}}_{1}(h)=I(f),}
wobei wir mit „
h
→
0
{\displaystyle h\to 0}
“ hier „
h
k
=
(
b
−
a
)
/
N
k
{\displaystyle h_{k}=(b-a)/N_{k}}
mit
N
k
∈
N
{\displaystyle N_{k}\in \mathbb {N} }
und
N
k
→
∞
{\displaystyle N_{k}\to \infty }
“ meinen. Die Entwicklung (8.23) soll nun numerisch dazu ausgenutzt werden, von einer endlichen Zahl berechneter Werte
J
1
(
h
k
)
,
k
=
0
,
1
,
…
,
n
{\displaystyle {\mathcal {J}}_{1}(h_{k}),k=0,1,\ldots ,n}
mit
0
<
h
n
<
h
n
−
1
<
…
<
h
0
{\displaystyle 0<h_{n}<h_{n-1}<\ldots <h_{0}}
auf einen noch genaueren Wert von
I
(
f
)
{\displaystyle I(f)}
als
J
1
(
h
n
)
{\displaystyle {\mathcal {J}}_{1}(h_{n})}
zu schließen.
Wir gehen dabei allgemeiner von einer beliebigen Funktion
T
(
h
)
{\displaystyle T(h)}
mit
h
>
0
{\displaystyle h>0}
aus, die mit gewissen Koeffizienten
α
i
∈
R
{\displaystyle \alpha _{i}\in \mathbb {R} }
(
i
=
0
,
1
,
…
,
r
)
{\displaystyle (i=0,1,\ldots ,r)}
und einer Zahl
γ
>
0
{\displaystyle \gamma >0}
die asymptotische Entwicklung der Ordnung
r
{\displaystyle r}
(8.24)
T
(
h
)
=
α
0
+
α
1
h
γ
+
α
2
h
2
γ
+
…
+
α
r
h
r
γ
+
O
(
h
(
r
+
1
)
γ
)
{\displaystyle T(h)=\alpha _{0}+\alpha _{1}h^{\gamma }+\alpha _{2}h^{2\gamma }+\ldots +\alpha _{r}h^{r\gamma }+{\mathcal {O}}(h^{(r+1)\gamma })}
für
h
→
0
{\displaystyle h\to 0}
besitzt und für die der Wert
lim
h
→
0
T
(
h
)
=
α
0
{\displaystyle \lim _{h\to 0}T(h)=\alpha _{0}}
gesucht ist. Typischerweise steht
T
(
h
)
{\displaystyle T(h)}
für ein numerisches Verfahren, das für einen gewählten Diskretisierungsparameter
h
>
0
{\displaystyle h>0}
einen Näherungswert für die gesuchte Größe
α
0
{\displaystyle \alpha _{0}}
liefert. Es sei also angenommen, dass
T
(
h
)
{\displaystyle T(h)}
zumindest für gewisse
h
>
0
{\displaystyle h>0}
berechnet werden kann, wie dies z. B. im Fall der Tangententrapezregel für
h
:=
(
b
−
a
)
/
N
{\displaystyle h:=(b-a)/N}
mit
N
∈
N
{\displaystyle N\in \mathbb {N} }
der Fall ist.
Wegen (8.24) hat man zunächst für
h
>
0
{\displaystyle h>0}
nur die Genauigkeit
T
(
h
)
−
α
0
=
O
(
h
γ
)
.
{\displaystyle T(h)-\alpha _{0}={\mathcal {O}}(h^{\gamma }).}
Es soll nun ein Verfahren vorgestellt werden, welches ohne großen Mehraufwand aus endlich vielen, bereits berechneten Werten
T
(
h
k
)
,
k
=
0
,
1
,
…
,
n
{\displaystyle T(h_{k}),k=0,1,\ldots ,n}
mit
0
<
h
n
<
h
n
−
1
<
…
<
h
0
{\displaystyle 0<h_{n}<h_{n-1}<\ldots <h_{0}}
einen genaueren Wert für die gesuchte Größe
α
0
{\displaystyle \alpha _{0}}
erzeugt. Setzt man
T
(
0
)
:=
α
0
{\displaystyle T(0):=\alpha _{0}}
, so extrapoliert dieses Verfahren also
T
{\displaystyle T}
auf den Wert
h
=
0
{\displaystyle h=0}
hin, so dass man auch von einem Extrapolationsverfahren spricht. Da die Koeffizienten
α
i
{\displaystyle \alpha _{i}}
in (8.24) oft nicht explizit bekannt sind oder nur unter einigem Aufwand zu berechnen sind, geht man dabei folgendermaßen vor:
man vernachlässigt den Restterm
O
(
h
(
r
+
1
)
γ
)
{\displaystyle {\mathcal {O}}(h^{(r+1)\gamma })}
in (8.24) und geht davon aus, dass sich
T
(
h
)
{\displaystyle T(h)}
ungefähr wie ein Polynom in
h
{\displaystyle h}
verhält,
man ersetzt das resultierende (i. A. nicht explizit bekannte) Polynom durch das Interpolationspolynom
P
0
,
…
,
n
∈
Π
n
{\displaystyle P_{0,\ldots ,n}\in \Pi _{n}}
zu den Stützpunkten
(
h
k
γ
,
T
(
h
k
)
)
,
k
=
0
,
1
,
…
,
n
{\displaystyle (h_{k}^{\gamma },T(h_{k})),k=0,1,\ldots ,n}
(schreibt man
T
(
h
γ
)
{\displaystyle T(h^{\gamma })}
statt
T
(
h
)
{\displaystyle T(h)}
, so sind dies mit
z
k
:=
h
k
γ
{\displaystyle z_{k}:=h_{k}^{\gamma }}
die Punkte
(
z
k
,
T
(
z
k
)
)
{\displaystyle (z_{k},T(z_{k}))}
) und
man verwendet den Wert
P
0
,
…
,
n
(
0
)
{\displaystyle P_{0,\ldots ,n}(0)}
als Näherung für den unbekannten Wert
α
0
{\displaystyle \alpha _{0}}
.
Im Zusammenhang mit der summierten Trapezregel wird diese Vorgehensweise als Romberg-Verfahren bezeichnet.
Wir gehen nun von der asymptotischen Entwicklung (8.24) von
T
(
h
)
{\displaystyle T(h)}
aus und es sei
P
0
,
…
,
n
∈
Π
n
{\displaystyle P_{0,\ldots ,n}\in \Pi _{n}}
das Interpolationspolynom zu den Stützpunkten
(8.25)
(
h
k
γ
,
T
(
h
k
)
)
k
=
0
,
1
,
…
,
n
.
{\displaystyle (h_{k}^{\gamma },T(h_{k}))\quad k=0,1,\ldots ,n.}
Da dieses nur an einer Stelle, der Stelle 0, ausgewertet werden soll, bietet sich das Neville-Schema zur Verwendung an, wobei hier
P
j
,
…
,
j
+
m
∈
Π
m
{\displaystyle P_{j,\ldots ,j+m}\in \Pi _{m}}
das Interpolationspolynom mit
(8.26)
P
j
,
…
,
j
+
m
(
h
k
γ
)
=
T
(
h
k
)
,
k
=
j
,
…
,
j
+
m
{\displaystyle P_{j,\ldots ,j+m}(h_{k}^{\gamma })=T(h_{k}),\quad k=j,\ldots ,j+m}
bezeichnet. Wir setzen dazu
(8.27)
T
j
,
…
,
j
+
m
:=
P
j
,
…
,
j
+
m
(
0
)
.
{\displaystyle T_{j,\ldots ,j+m}:=P_{j,\ldots ,j+m}(0).}
Satz 6.5 liefert damit
(8.28)
T
j
=
P
j
(
0
)
=
T
(
h
j
)
{\displaystyle T_{j}=P_{j}(0)=T(h_{j})}
sowie
T
j
,
…
,
j
+
m
=
−
h
j
γ
T
j
+
1
,
…
,
j
+
m
+
h
j
+
m
γ
T
j
,
…
,
j
+
m
−
1
h
j
+
m
γ
−
h
j
γ
=
T
j
+
1
,
…
,
j
+
m
−
h
j
+
m
γ
T
j
+
1
,
…
,
j
+
m
−
T
j
,
…
,
j
+
m
−
1
h
j
+
m
γ
−
h
j
γ
{\displaystyle T_{j,\ldots ,j+m}={\frac {-h_{j}^{\gamma }T_{j+1,\ldots ,j+m}+h_{j+m}^{\gamma }T_{j,\ldots ,j+m-1}}{h_{j+m}^{\gamma }-h_{j}^{\gamma }}}=T_{j+1,\ldots ,j+m}-h_{j+m}^{\gamma }{\frac {T_{j+1,\ldots ,j+m}-T_{j,\ldots ,j+m-1}}{h_{j+m}^{\gamma }-h_{j}^{\gamma }}}}
(8.29)
=
T
j
+
1
,
…
,
j
+
m
+
T
j
+
1
,
…
,
j
+
m
−
T
j
,
…
,
j
+
m
−
1
(
h
j
h
j
+
m
)
γ
−
1
(
j
,
m
∈
N
0
)
.
{\displaystyle =T_{j+1,\ldots ,j+m}+{\frac {T_{j+1,\ldots ,j+m}-T_{j,\ldots ,j+m-1}}{\left({\frac {h_{j}}{h_{j+m}}}\right)^{\gamma }-1}}\quad (j,m\in \mathbb {N} _{0}).}
Das Schema von Neville geht damit in das folgende Extrapolationstableau über, welches zeilenweise aufgebaut wird:
T
0
=
T
(
h
0
)
↘
T
1
=
T
(
h
1
)
→
T
01
↘
↘
T
2
=
T
(
h
2
)
→
T
12
→
T
012
↘
↘
↘
T
3
=
T
(
h
3
)
→
T
23
→
T
123
→
T
0123
⋮
⋮
⋮
⋮
⋱
{\displaystyle {\begin{matrix}T_{0}=T(h_{0})&&&&&&&\\&\searrow &&&&&&\\T_{1}=T(h_{1})&\to &T_{01}&&&&&\\&\searrow &&\searrow &&&&\\T_{2}=T(h_{2})&\to &T_{12}&\to &T_{012}&&&\\&\searrow &&\searrow &&\searrow &&\\T_{3}=T(h_{3})&\to &T_{23}&\to &T_{123}&\to &T_{0123}&\\\vdots &&\vdots &&\vdots &&\vdots &\ddots \end{matrix}}}
Für die summierte Trapezregel
J
1
(
h
)
{\displaystyle {\mathcal {J}}_{1}(h)}
gilt gemäß Satz 8.23 eine Entwicklung der Form (8.24) mit
γ
=
2
{\displaystyle \gamma =2}
. Für die Schrittweiten
h
0
:=
b
−
a
,
h
1
:=
h
0
2
=
b
−
a
2
{\displaystyle h_{0}:=b-a,\quad h_{1}:={\frac {h_{0}}{2}}={\frac {b-a}{2}}}
erhält man die Werte
T
0
=
J
1
(
h
0
)
=
b
−
a
2
[
f
(
a
)
+
f
(
b
)
]
,
{\displaystyle T_{0}={\mathcal {J}}_{1}(h_{0})={\frac {b-a}{2}}[f(a)+f(b)],}
T
1
=
J
1
(
h
0
)
=
b
−
a
2
[
f
(
a
)
+
2
f
(
a
+
b
2
)
+
f
(
b
)
]
{\displaystyle T_{1}={\mathcal {J}}_{1}(h_{0})={\frac {b-a}{2}}\left[f(a)+2f\left({\frac {a+b}{2}}\right)+f(b)\right]}
und damit
T
01
=
T
1
+
T
1
−
T
0
(
h
0
h
1
)
2
−
1
=
b
−
a
4
[
f
(
a
)
+
2
f
(
a
+
b
2
)
+
f
(
b
)
]
+
b
−
a
4
⋅
2
f
(
a
+
b
2
)
−
f
(
a
)
−
f
(
b
)
4
−
1
{\displaystyle T_{01}=T_{1}+{\frac {T_{1}-T_{0}}{\left({\frac {h_{0}}{h_{1}}}\right)^{2}-1}}={\frac {b-a}{4}}\left[f(a)+2f\left({\frac {a+b}{2}}\right)+f(b)\right]+{\frac {b-a}{4}}\cdot {\frac {2f\left({\frac {a+b}{2}}\right)-f(a)-f(b)}{4-1}}}
=
b
−
a
4
[
2
3
f
(
a
)
+
8
3
f
(
a
+
b
2
)
+
2
3
f
(
b
)
]
=
b
−
a
6
[
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
]
.
{\displaystyle ={\frac {b-a}{4}}\left[{\frac {2}{3}}f(a)+{\frac {8}{3}}f\left({\frac {a+b}{2}}\right)+{\frac {2}{3}}f(b)\right]={\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right].}
Der aus den beiden Auswertungen
T
0
{\displaystyle T_{0}}
und
T
1
{\displaystyle T_{1}}
der summierten Trapezregel ermittelte Wert
T
01
{\displaystyle T_{01}}
entspricht somit dem der Simpson-Regel für
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
.
Im folgenden Satz wird die Größenordnung des Fehlers
T
j
,
…
,
j
+
m
−
α
0
{\displaystyle T_{j,\ldots ,j+m}-\alpha _{0}}
angegeben. Diese Fehlerbetrachtung macht deutlich, dass sich die Anwendung des hier untersuchten Extrapolationsverfahrens lohnt. Als Hilfsmittel verwenden wir das nachstehende Lemma.
Es seien
L
k
∈
Π
m
{\displaystyle L_{k}\in \Pi _{m}}
(
k
=
0
,
1
,
…
,
m
)
{\displaystyle (k=0,1,\ldots ,m)}
die Lagrangeschen Basispolynome zu Stützstellen
x
k
{\displaystyle x_{k}}
(
k
=
0
,
1
,
…
,
m
)
{\displaystyle (k=0,1,\ldots ,m)}
mit
x
k
≠
x
i
{\displaystyle x_{k}\neq x_{i}}
(
k
≠
i
)
{\displaystyle (k\neq i)}
. Dann gilt
(8.30)
∑
k
=
0
m
L
k
(
0
)
x
k
j
=
{
1
f
u
¨
r
j
=
0
,
0
f
u
¨
r
1
≤
j
≤
m
,
(
−
1
)
m
x
0
x
1
⋯
x
m
f
u
¨
r
j
=
m
+
1.
{\displaystyle \sum _{k=0}^{m}L_{k}(0)x_{k}^{j}={\begin{cases}1&f{\ddot {u}}r\ j=0,\\0&f{\ddot {u}}r\ 1\leq j\leq m,\\(-1)^{m}x_{0}x_{1}\cdots x_{m}&f{\ddot {u}}r\ j=m+1.\end{cases}}}
Für
0
≤
j
≤
m
{\displaystyle 0\leq j\leq m}
ist offenbar
p
(
x
)
:=
x
j
{\displaystyle p(x):=x^{j}}
das Interpolationspolynom zu den Punkten
(
x
k
,
x
k
j
)
,
k
=
0
,
1
,
…
,
m
{\displaystyle (x_{k},x_{k}^{j}),k=0,1,\ldots ,m}
und daher gemäß (6.7)
p
(
x
)
=
x
j
=
∑
k
=
0
m
x
k
j
L
k
(
x
)
.
{\displaystyle p(x)=x^{j}=\sum _{k=0}^{m}x_{k}^{j}L_{k}(x).}
Setzen wir
x
=
0
{\displaystyle x=0}
, so folgt die Behauptung für die ersten beiden Fälle in (8.30). Für den Fall
j
=
m
+
1
{\displaystyle j=m+1}
betrachten wir das Polynom
q
(
x
)
:=
x
m
+
1
−
∑
k
=
0
m
x
k
m
+
1
L
k
(
x
)
,
{\displaystyle q(x):=x^{m+1}-\sum _{k=0}^{m}x_{k}^{m+1}L_{k}(x),}
welches wegen
L
k
∈
Π
m
{\displaystyle L_{k}\in \Pi _{m}}
den Grad
m
+
1
{\displaystyle m+1}
, den führenden Koeffizienten 1 und die Nullstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
m
)
{\displaystyle (i=0,1,\ldots ,m)}
hat, so dass insbesondere
q
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
m
)
{\displaystyle q(x)=(x-x_{0})(x-x_{1})\cdots (x-x_{m})}
gilt. Speziell hat man somit
∑
k
=
0
m
L
k
(
0
)
x
k
m
+
1
=
−
q
(
0
)
=
(
−
1
)
m
x
0
x
1
⋯
x
m
.
{\displaystyle \sum _{k=0}^{m}L_{k}(0)x_{k}^{m+1}=-q(0)=(-1)^{m}x_{0}x_{1}\cdots x_{m}.}
Es sei
T
(
h
)
{\displaystyle T(h)}
mit
h
>
0
{\displaystyle h>0}
eine Funktion mit der asymptotischen Entwicklung (8.24) für ein
γ
>
0
{\displaystyle \gamma >0}
und
r
∈
N
{\displaystyle r\in \mathbb {N} }
. Weiter sei
(
h
k
)
{\displaystyle (h_{k})}
eine Folge von Schrittweiten, so dass mit einer Startschrittweite
h
0
>
0
{\displaystyle h_{0}>0}
gilt:
(8.31)
h
k
:=
h
0
/
n
k
(
k
∈
N
0
)
{\displaystyle h_{k}:=h_{0}/n_{k}\quad (k\in \mathbb {N} _{0})}
mit
1
=
n
0
≤
n
1
≤
n
2
≤
…
.
{\displaystyle 1=n_{0}\leq n_{1}\leq n_{2}\leq \ldots .}
Schließlich sei
P
j
,
…
,
j
+
m
∈
Π
m
{\displaystyle P_{j,\ldots ,j+m}\in \Pi _{m}}
das Interpolationspolynom mit (8.26) und
T
j
,
…
,
j
+
m
{\displaystyle T_{j,\ldots ,j+m}}
wie in (8.27). Dann genügt der Fehler
T
j
,
…
,
j
+
m
−
α
0
{\displaystyle T_{j,\ldots ,j+m}-\alpha _{0}}
für
0
≤
m
≤
r
−
1
{\displaystyle 0\leq m\leq r-1}
der asymptotischen Entwicklung
(8.32)
T
j
,
…
,
j
+
m
−
α
0
=
(
−
1
)
m
α
m
+
1
n
j
γ
⋯
n
j
+
m
γ
h
0
(
m
+
1
)
γ
+
O
(
h
0
(
m
+
2
)
γ
)
.
{\displaystyle T_{j,\ldots ,j+m}-\alpha _{0}=(-1)^{m}{\frac {\alpha _{m+1}}{n_{j}^{\gamma }\cdots n_{j+m}^{\gamma }}}h_{0}^{(m+1)\gamma }+{\mathcal {O}}(h_{0}^{(m+2)\gamma }).}
Da sich die Indizes in (8.32) auf eine Numerierung der Stützpunkte beziehen und wir den
j
{\displaystyle j}
-ten als 0-ten bezeichnen können, können wir o. B. d. A.
j
=
0
{\displaystyle j=0}
annehmen. Gemäß der Lagrangeschen Darstellung des Interpolationspolynoms
P
0
,
…
,
m
{\displaystyle P_{0,\ldots ,m}}
gilt dann
P
0
,
…
,
m
(
h
γ
)
=
∑
k
=
0
m
T
(
h
k
)
L
k
(
h
γ
)
=
∑
k
=
0
m
T
(
h
k
)
[
∏
j
=
0
j
≠
k
m
h
−
h
j
γ
h
k
γ
−
h
j
γ
]
,
h
∈
R
{\displaystyle P_{0,\ldots ,m}(h^{\gamma })=\sum _{k=0}^{m}T(h_{k})L_{k}(h^{\gamma })=\sum _{k=0}^{m}T(h_{k})\left[\prod _{j=0 \atop j\neq k}^{m}{\frac {h-h_{j}^{\gamma }}{h_{k}^{\gamma }-h_{j}^{\gamma }}}\right],\quad h\in \mathbb {R} }
und somit
(8.33)
T
0
,
…
,
m
=
P
0
,
…
,
m
(
0
)
=
∑
k
=
0
m
c
m
,
k
T
(
h
k
)
{\displaystyle T_{0,\ldots ,m}=P_{0,\ldots ,m}(0)=\sum _{k=0}^{m}c_{m,k}T(h_{k})}
für
c
m
,
k
:=
L
k
(
0
)
=
∏
j
=
0
j
≠
k
m
h
j
γ
h
j
γ
−
h
k
γ
.
{\displaystyle c_{m,k}:=L_{k}(0)=\prod _{j=0 \atop j\neq k}^{m}{\frac {h_{j}^{\gamma }}{h_{j}^{\gamma }-h_{k}^{\gamma }}}.}
Nun folgt wegen
m
≤
r
−
1
{\displaystyle m\leq r-1}
aus (8.24)
(8.34)
T
(
h
k
)
=
∑
s
=
0
m
+
1
α
s
h
k
s
γ
+
O
(
h
k
(
m
+
2
)
γ
)
.
{\displaystyle T(h_{k})=\sum _{s=0}^{m+1}\alpha _{s}h_{k}^{s\gamma }+{\mathcal {O}}(h_{k}^{(m+2)\gamma }).}
Des Weiteren schließt man mit Lemma 8.25
(8.35)
∑
k
=
0
m
c
m
,
k
h
k
s
γ
=
∑
k
=
0
m
L
k
(
0
)
h
k
s
γ
=
{
1
für
s
=
0
,
0
für
1
≤
s
≤
m
,
(
−
1
)
m
h
0
γ
h
1
γ
⋯
h
m
γ
für
s
=
m
+
1.
{\displaystyle \sum _{k=0}^{m}c_{m,k}h_{k}^{s\gamma }=\sum _{k=0}^{m}L_{k}(0)h_{k}^{s\gamma }={\begin{cases}1&{\mbox{für }}s=0,\\0&{\mbox{für }}1\leq s\leq m,\\(-1)^{m}h_{0}^{\gamma }h_{1}^{\gamma }\cdots h_{m}^{\gamma }&{\mbox{für }}s=m+1.\end{cases}}}
Setzt man die beiden Beziehungen (8.34) und (8.35) in (8.33) ein, so bekommt man schließlich, da die
c
m
,
k
{\displaystyle c_{m,k}}
von
h
{\displaystyle h}
unabhängig sind,
T
0
,
…
,
m
=
∑
k
=
0
m
c
m
,
k
[
∑
s
=
0
m
+
1
α
s
h
k
s
γ
+
O
(
h
k
(
m
+
2
)
γ
)
]
=
∑
s
=
0
m
+
1
α
s
[
∑
k
=
0
m
c
m
,
k
h
k
s
γ
]
+
∑
k
=
0
m
c
m
,
k
O
(
h
k
(
m
+
2
)
γ
)
{\displaystyle T_{0,\ldots ,m}=\sum _{k=0}^{m}c_{m,k}\left[\sum _{s=0}^{m+1}\alpha _{s}h_{k}^{s\gamma }+{\mathcal {O}}(h_{k}^{(m+2)\gamma })\right]=\sum _{s=0}^{m+1}\alpha _{s}\left[\sum _{k=0}^{m}c_{m,k}h_{k}^{s\gamma }\right]+\sum _{k=0}^{m}c_{m,k}{\mathcal {O}}(h_{k}^{(m+2)\gamma })}
=
α
0
+
(
−
1
)
m
α
m
+
1
h
0
γ
h
1
γ
⋯
h
m
γ
+
O
(
h
0
(
m
+
2
)
γ
)
.
{\displaystyle =\alpha _{0}+(-1)^{m}\alpha _{m+1}h_{0}^{\gamma }h_{1}^{\gamma }\cdots h_{m}^{\gamma }+{\mathcal {O}}(h_{0}^{(m+2)\gamma }).}
q.e.d.
Der Satz besagt, dass man beim Übergang von
m
{\displaystyle m}
zu
m
+
1
{\displaystyle m+1}
, d. h. bei Erhöhung der Spaltenzahl in dem Extrapolationstableau um 1, im Prinzip die Ordnung
γ
{\displaystyle \gamma }
gewinnt. Diese Sichtweise ist allerdings zu optimistisch, da die Restterme der asymptotischen Entwicklung, die sich hinter
O
(
h
0
(
m
+
2
)
γ
)
{\displaystyle {\mathcal {O}}(h_{0}^{(m+2)\gamma })}
verbergen, nicht bekannt sind und groß werden können.
Es bietet sich also der folgende Algorithmus an:
(0) Wähle
h
0
>
0
{\displaystyle h_{0}>0}
, eine Folge
h
k
,
k
=
1
,
2
,
…
{\displaystyle h_{k},k=1,2,\ldots }
wie in (8.31) und ein
ε
>
0
{\displaystyle \varepsilon >0}
. Setze
j
:=
0
{\displaystyle j:=0}
.
(1) Berechne
T
j
:=
T
(
h
j
)
{\displaystyle T_{j}:=T(h_{j})}
.
(2) Berechne
T
k
,
…
,
j
{\displaystyle T_{k,\ldots ,j}}
für
k
=
j
−
1
,
j
−
2
,
…
,
0
{\displaystyle k=j-1,j-2,\ldots ,0}
nach der Formel
T
k
,
…
,
j
=
T
k
+
1
,
…
,
j
+
T
k
+
1
,
…
,
j
−
T
k
,
…
,
j
−
1
(
h
k
h
j
)
γ
−
1
{\displaystyle T_{k,\ldots ,j}=T_{k+1,\ldots ,j}+{\frac {T_{k+1,\ldots ,j}-T_{k,\ldots ,j-1}}{\left({\frac {h_{k}}{h_{j}}}\right)^{\gamma }-1}}}
(3) Falls „der Aufwand zu groß wird“ oder
|
T
0
,
…
,
j
−
T
0
,
…
,
j
−
1
|
|
T
0
,
…
,
j
|
{\displaystyle {\frac {|T_{0,\ldots ,j}-T_{0,\ldots ,j-1}|}{|T_{0,\ldots ,j}|}}}
gilt, breche ab. (
T
0
,
…
,
j
{\displaystyle T_{0,\ldots ,j}}
ist Näherungswert für
α
0
{\displaystyle \alpha _{0}}
.)
(4) Setze
j
:=
j
+
1
{\displaystyle j:=j+1}
und gehe nach (1).
Man bricht das Extrapolationsverfahren also ab, wenn der Aufwand zur Erzeugung einer neuen Zeile im Extrapolationsschema, den man meistens, wie z. B. für das summierte Trapezverfahren, genau angeben kann, zu groß wird oder die relative Abweichung zweier aufeinanderfolgender Diagonalelemente klein genug wird. In der Praxis ist es jedoch auch möglich, dass aufgrund von Rundungsfehlern Divergenz eintritt, so dass auf früher berechnete Werte im Schema zurückgegriffen werden muss.
Häufig angewandte Schrittweitenfolgen
(
h
k
)
{\displaystyle (h_{k})}
für (8.31) in diesem Zusammenhang sind die Romberg-Folge
(8.36)
n
k
:=
2
k
,
h
k
=
h
k
−
1
2
=
h
0
2
k
(
k
∈
N
)
,
{\displaystyle n_{k}:=2^{k},\quad h_{k}={\frac {h_{k-1}}{2}}={\frac {h_{0}}{2^{k}}}\quad (k\in \mathbb {N} ),}
die durch
n
1
:=
2
,
n
2
:=
3
,
n
3
:=
4
,
n
j
:=
2
n
j
−
2
(
j
≥
4
)
{\displaystyle n_{1}:=2,\quad n_{2}:=3,\quad n_{3}:=4,\quad n_{j}:=2n_{j-2}\quad (j\geq 4)}
definierte Bulirsch-Folge
h
1
=
h
0
2
,
h
2
=
h
0
3
,
h
3
=
h
0
4
,
h
4
=
h
0
6
,
h
5
=
h
0
8
,
h
6
=
h
0
12
,
h
7
=
h
0
16
,
…
{\displaystyle h_{1}={\frac {h_{0}}{2}},\quad h_{2}={\frac {h_{0}}{3}},\quad h_{3}={\frac {h_{0}}{4}},\quad h_{4}={\frac {h_{0}}{6}},\quad h_{5}={\frac {h_{0}}{8}},\quad h_{6}={\frac {h_{0}}{12}},\quad h_{7}={\frac {h_{0}}{16}},\quad \ldots }
und die harmonische Folge
n
k
−
1
:=
k
(
k
=
1
,
2
,
…
)
,
h
k
=
h
0
k
+
1
(
k
∈
N
)
.
{\displaystyle n_{k-1}:=k\quad (k=1,2,\ldots ),\quad h_{k}={\frac {h_{0}}{k+1}}\quad (k\in \mathbb {N} ).}
Insbesondere erhält man für die Romberg-Folge
(
j
:=
0
)
{\displaystyle (j:=0)}
:
Unter den Voraussetzungen von Satz 8.26 gilt für die Romberg-Folge (8.36) mit
h
0
:=
(
b
−
a
)
/
N
{\displaystyle h_{0}:=(b-a)/N}
und
0
≤
m
≤
r
−
1
{\displaystyle 0\leq m\leq r-1}
T
0
,
…
,
m
−
α
0
=
[
(
−
1
)
m
2
γ
m
(
m
+
1
)
/
2
α
m
+
1
]
h
0
(
m
+
1
)
γ
+
O
(
h
0
(
m
+
2
)
γ
)
.
{\displaystyle T_{0,\ldots ,m}-\alpha _{0}=\left[{\frac {(-1)^{m}}{2^{\gamma m(m+1)/2}}}\alpha _{m+1}\right]h_{0}^{(m+1)\gamma }+{\mathcal {O}}(h_{0}^{(m+2)\gamma }).}
Im Fall (8.36) hat man mit
n
0
=
1
{\displaystyle n_{0}=1}
n
0
γ
⋯
n
m
γ
=
2
0
γ
⋅
2
1
γ
⋅
2
2
γ
⋯
2
m
γ
=
2
γ
∑
k
=
0
m
k
=
2
γ
m
(
m
+
1
)
/
2
,
{\displaystyle n_{0}^{\gamma }\cdots n_{m}^{\gamma }=2^{0\gamma }\cdot 2^{1\gamma }\cdot 2^{2\gamma }\cdots 2^{m\gamma }=2^{\gamma \sum _{k=0}^{m}k}=2^{\gamma m(m+1)/2},}
so dass die Behauptung unmittelbar aus Satz 8.26 folgt.
q.e.d.
Wir betrachten nun nochmals die summierte Trapezregel als Beispiel.
(1) Korollar 8.27 wollen wir auf die summierte Trapezregel mit
m
:=
2
{\displaystyle m:=2}
(und wegen der Forderung
m
≤
r
−
1
{\displaystyle m\leq r-1}
für
r
=
3
{\displaystyle r=3}
) anwenden. Nach Satz 8.23 ist dann
γ
=
2
{\displaystyle \gamma =2}
. Weiter sei
f
∈
C
8
[
a
,
b
]
{\displaystyle f\in C^{8}[a,b]}
vorausgesetzt. Korollar 8.27 liefert mit diesen Setzungen
(8.37)
T
012
−
I
(
f
)
=
[
(
−
1
)
2
2
2
⋅
2
(
2
+
1
)
/
2
α
2
+
1
]
h
0
(
2
+
1
)
2
+
O
(
h
0
(
m
+
2
)
2
)
=
α
3
64
h
0
6
+
O
(
h
0
8
)
,
{\displaystyle T_{012}-{\mathcal {I}}(f)=\left[{\frac {(-1)^{2}}{2^{2\cdot 2(2+1)/2}}}\alpha _{2+1}\right]h_{0}^{(2+1)2}+{\mathcal {O}}(h_{0}^{(m+2)2})={\frac {\alpha _{3}}{64}}h_{0}^{6}+{\mathcal {O}}(h_{0}^{8}),}
wobei
T
012
{\displaystyle T_{012}}
mit dem Neville-Schema
T
0
=
J
1
(
h
0
)
↘
T
1
=
J
1
(
h
0
/
2
)
→
T
01
↘
↘
T
2
=
J
1
(
h
0
/
4
)
→
T
12
→
T
012
{\displaystyle {\begin{matrix}T_{0}={\mathcal {J}}_{1}(h_{0})&&&&\\&\searrow &&&\\T_{1}={\mathcal {J}}_{1}(h_{0}/2)&\rightarrow &T_{01}&&\\&\searrow &&\searrow &\\T_{2}={\mathcal {J}}_{1}(h_{0}/4)&\rightarrow &T_{12}&\rightarrow &T_{012}\end{matrix}}}
berechnet wird. Man beachte dabei, dass man bei der Berechnung von
T
j
+
1
{\displaystyle T_{j+1}}
(
j
∈
N
0
)
{\displaystyle (j\in \mathbb {N} _{0})}
den zuvor ermittelten Wert
T
j
{\displaystyle T_{j}}
verwenden kann und nur zusätzlich Funktionsauswertungen für die Mittelpunkte der sich aus der zu
T
j
{\displaystyle T_{j}}
gehörenden Zerlegung von
[
a
,
b
]
{\displaystyle [a,b]}
ergebenden Intervalle benötigt. Somit verlangt die Berechnung von
T
0
,
T
1
{\displaystyle T_{0},T_{1}}
und
T
2
{\displaystyle T_{2}}
genauso viele Funktionsauswertungen wie die direkte Berechnung von
T
2
{\displaystyle T_{2}}
. Für letzteren Wert alleine hat man aber im Vergleich zu (8.37) gemäß Satz 8.20 für
f
∈
C
2
[
a
,
b
]
{\displaystyle f\in C^{2}[a,b]}
mit einem
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
einen Fehler der Größe
O
(
h
0
2
)
{\displaystyle {\mathcal {O}}(h_{0}^{2})}
:
T
2
−
I
(
f
)
=
J
1
(
h
0
/
4
)
−
I
(
f
)
=
b
−
a
12
⋅
4
2
h
0
2
f
″
(
ξ
)
=
(
b
−
a
)
f
″
(
ξ
)
/
3
64
h
0
2
.
{\displaystyle T_{2}-{\mathcal {I}}(f)={\mathcal {J}}_{1}(h_{0}/4)-{\mathcal {I}}(f)={\frac {b-a}{12\cdot 4^{2}}}h_{0}^{2}f''(\xi )={\frac {(b-a)f''(\xi )/3}{64}}h_{0}^{2}.}
(2) (Bader) Es soll das Integral
∫
0
π
/
2
5
e
2
x
cos
(
x
)
e
π
−
2
d
x
=
1
{\displaystyle \int \limits _{0}^{\pi /2}5{\frac {e^{2x}\cos(x)}{e^{\pi }-2}}\,dx=1}
näherungsweise mit der summierten Trapezregel und dem Extrapolationsverfahren mit der Romberg-Folge (8.36) und
h
0
:=
π
/
2
{\displaystyle h_{0}:=\pi /2}
berechnet werden. Es ergibt sich bei 12-stelliger Rechnung das folgende Extrapolationstableau:
0.185
755
068
924
0.724
727
335
089
0.904
384
757
145
0.925
565
035
158
0.992
510
935
182
0.999
386
013
717
0.981
021
630
069
0.999
507
161
706
0.999
973
576
808
0.999
998
776
222
0.995
232
017
388
0.999
968
813
161
0.999
999
589
925
1.000
000
002
83
0.999
806
537
974
0.999
998
044
836
0.999
999
993
614
1.000
000
000
02
{\displaystyle {\begin{matrix}0.185\,755\,068\,924&&&\\0.724\,727\,335\,089&0.904\,384\,757\,145&&\\0.925\,565\,035\,158&0.992\,510\,935\,182&0.999\,386\,013\,717&\\0.981\,021\,630\,069&0.999\,507\,161\,706&0.999\,973\,576\,808&0.999\,998\,776\,222\\0.995\,232\,017\,388&0.999\,968\,813\,161&0.999\,999\,589\,925&1.000\,000\,002\,83\\0.999\,806\,537\,974&0.999\,998\,044\,836&0.999\,999\,993\,614&1.000\,000\,000\,02\end{matrix}}}
Der (in der Tabelle nicht mehr einfügbare) Wert des Diagonalelementes in der 5. Spalte beträgt
1.000
000
084
6
{\displaystyle 1.000\,000\,084\,6}
. Er ist offenbar ungenauer als die beiden untersten Werte in der 4. Spalte, wobei für den untersten allerdings auch die summierte Trapezregel einmal mehr ausgewertet werden musste. (Man kann auch zeigen, was für die erste Spalte z. B. aus Satz 8.20 folgt, dass die Werte in jeder einzelnen Spalte des Extrapolationsschemas, d. h. für konstantes
m
{\displaystyle m}
und
j
→
∞
{\displaystyle j\to \infty }
, gegen den gesuchten Wert
α
0
{\displaystyle \alpha _{0}}
konvergieren.) Für eine Diskussion über ein geeignetes Abbruchkriterium verweisen wir auf Deuflhard/Hohmann.
In diesem Abschnitt betrachten wir Quadraturformeln für gewichtete Integrale des Typs
(8.38)
I
(
f
)
:=
∫
a
b
w
(
x
)
f
(
x
)
d
x
,
f
∈
C
[
a
,
b
]
,
{\displaystyle {\mathcal {I}}(f):=\int \limits _{a}^{b}w(x)f(x)\,dx,\quad f\in C[a,b],}
wobei das Intervall
[
a
,
b
]
{\displaystyle [a,b]}
hier halbunendlich oder unendlich, d. h.
a
:=
−
∞
{\displaystyle a:=-\infty }
und/oder
b
:=
∞
{\displaystyle b:=\infty }
sein darf und
w
:
(
a
,
b
)
→
R
{\displaystyle w:(a,b)\to \mathbb {R} }
eine Gewichtsfunktion mit den folgenden Eigenschaften ist:
w
(
x
)
>
0
,
x
∈
(
a
,
b
)
,
{\displaystyle w(x)>0,\quad x\in (a,b),}
es existieren die Momente
μ
k
:=
∫
a
b
w
(
x
)
x
k
d
x
,
k
=
0
,
1
,
2
,
…
.
{\displaystyle \mu _{k}:=\int \limits _{a}^{b}w(x)x^{k}\,dx,\quad k=0,1,2,\ldots .}
Häufig in diesem Zusammenhang auftretende Gewichtsfunktionen sind in der folgenden Tabelle wiedergegeben, wobei auch der zuvor betrachtete Fall
w
≡
1
{\displaystyle w\equiv 1}
von Interesse ist:
Fehler beim Parsen (Konvertierungsfehler. Der Server („https://wikimedia.org/api/rest_“) hat berichtet: „Cannot get mml. TeX parse error: Bracket argument to \\ must be a dimension“): {\displaystyle {\begin{array}{|c|c|}\hline {\text{Intervall}}&{\text{Gewichtsfunktion }}w(x)\\\hline [-1,1]&1\\[-1,1]&1/{\sqrt {1-x^{2}}}\\[0,\infty )&e^{-x}\\(-\infty ,\infty )&e^{-x^{2}}\\[0,\infty )&e^{-x^{2}}x^{\alpha },\ \alpha >-1\\\hline \end{array}}}
Wir definieren in diesem Zusammenhang auf dem Raum aller Polynome
Π
{\displaystyle \Pi }
das durch
w
{\displaystyle w}
induzierte Skalarprodukt
(8.39)
⟨
f
,
g
⟩
:=
∫
a
b
w
(
x
)
f
(
x
)
g
(
x
)
d
x
,
f
,
g
∈
Π
.
{\displaystyle \langle f,g\rangle :=\int \limits _{a}^{b}w(x)f(x)g(x)\,dx,\quad f,g\in \Pi .}
Das Integral in (8.39) existiert offenbar unter den Voraussetzungen an
w
{\displaystyle w}
. Für alle
f
,
g
,
h
∈
Π
{\displaystyle f,g,h\in \Pi }
gilt weiter (man verifiziere dies)
⟨
f
,
f
⟩
≥
0
,
⟨
f
,
f
⟩
=
0
⇔
f
≡
0
,
⟨
f
,
g
⟩
=
⟨
g
,
f
⟩
,
{\displaystyle \langle f,f\rangle \geq 0,\quad \langle f,f\rangle =0\Leftrightarrow f\equiv 0,\quad \langle f,g\rangle =\langle g,f\rangle ,}
⟨
α
f
+
β
g
,
h
⟩
=
α
⟨
f
,
h
⟩
+
β
⟨
g
,
h
⟩
=
⟨
h
,
α
f
+
β
g
⟩
,
α
,
β
∈
R
.
{\displaystyle \langle \alpha f+\beta g,h\rangle =\alpha \langle f,h\rangle +\beta \langle g,h\rangle =\langle h,\alpha f+\beta g\rangle ,\quad \alpha ,\beta \in \mathbb {R} .}
Insbesondere ist also die Abbildung
⟨
⋅
,
⋅
⟩
:
Π
×
Π
→
R
{\displaystyle \langle \cdot ,\cdot \rangle :\Pi \times \Pi \to \mathbb {R} }
in beiden Eingängen linear. Wir verwenden ferner die durch das Skalarprodukt
⟨
⋅
,
⋅
⟩
{\displaystyle \langle \cdot ,\cdot \rangle }
induzierte Norm auf
Π
{\displaystyle \Pi }
(8.40)
‖
f
‖
:=
⟨
f
,
f
⟩
1
/
2
=
{
∫
a
b
w
(
x
)
f
2
(
x
)
d
x
}
1
/
2
,
f
∈
Π
.
{\displaystyle \|f\|:=\langle f,f\rangle ^{1/2}=\left\{\int \limits _{a}^{b}w(x)f^{2}(x)dx\right\}^{1/2},\quad f\in \Pi .}
Ziel ist es nun wieder, zur numerischen Berechnung von
I
(
f
)
{\displaystyle {\mathcal {I}}(f)}
eine Quadraturformel
(8.41)
I
n
(
f
)
:=
∑
i
=
0
n
σ
i
f
(
x
i
)
{\displaystyle {\mathcal {I}}_{n}(f):=\sum _{i=0}^{n}\sigma _{i}f(x_{i})}
herzuleiten. (Man beachte, dass hier der Faktor
b
−
a
{\displaystyle b-a}
vor der Summe fehlt.) Und zwar soll eine interpolatorische Quadraturformel entwickelt werden, für welche bei geeigneter Wahl der Stützstellen
x
i
{\displaystyle x_{i}}
und der Gewichte
σ
i
{\displaystyle \sigma _{i}}
der Genauigkeitsgrad möglichst hoch ist, welche also Polynome bis zu einem möglichst hohen Grad exakt integriert. Man betrachte dazu die Aussagen in Satz 8.15 über den Quadraturfehler. Die Begriffe Genauigkeitsgrad und interpolatorische Quadraturformel sind hierbei analog zu den Definitionen 8.2 und 8.8 auf Integrale mit Gewichten zu übertragen.
Zunächst einmal stellen wir fest, dass man in (8.41) insgesamt 2n+2 Parameter
x
i
{\displaystyle x_{i}}
und
σ
i
{\displaystyle \sigma _{i}}
zur Verfügung hat, was der Anzahl der Koeffizienten eines Polynoms vom Grad
2
n
+
1
{\displaystyle 2n+1}
entspricht. In der Tat werden wir zeigen, dass eine Quadraturformel mit diesem Genauigkeitsgrad existiert. Quadraturformeln mit einem höheren Genauigkeitsgrad kann es nicht geben. Denn wäre
I
n
{\displaystyle {\mathcal {I}}_{n}}
eine Quadraturformel mit Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und Gewichten
σ
i
{\displaystyle \sigma _{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und hätte
I
n
{\displaystyle {\mathcal {I}}_{n}}
den Genauigkeitsgrad
2
n
+
2
{\displaystyle 2n+2}
, so folgte insbesondere für das Polynom
p
(
x
)
:=
∏
i
=
0
n
(
x
−
x
i
)
2
∈
Π
2
n
+
2
{\displaystyle p(x):=\prod _{i=0}^{n}(x-x_{i})^{2}\in \Pi _{2n+2}}
p
(
x
i
)
=
0
{\displaystyle p(x_{i})=0}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und daher
I
n
(
p
)
=
0
=
I
(
p
)
{\displaystyle {\mathcal {I}}_{n}(p)=0={\mathcal {I}}(p)}
. Wegen
p
(
x
)
>
0
,
x
∈
(
x
i
,
x
i
+
1
)
,
i
=
0
,
1
,
…
,
n
−
1
{\displaystyle p(x)>0,\quad x\in (x_{i},x_{i+1}),\quad i=0,1,\ldots ,n-1}
ist jedoch
I
(
p
)
>
0
{\displaystyle {\mathcal {I}}(p)>0}
. Wir können weiter schließen:
Ist
I
n
{\displaystyle {\mathcal {I}}_{n}}
eine Quadraturformel mit Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und Gewichten
σ
i
{\displaystyle \sigma _{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und hat
I
n
{\displaystyle {\mathcal {I}}_{n}}
den Genauigkeitsgrad
2
n
+
1
{\displaystyle 2n+1}
, so gilt
⟨
p
,
p
n
+
1
⟩
=
0
,
p
∈
Π
n
{\displaystyle \langle p,p_{n+1}\rangle =0,\quad p\in \Pi _{n}}
für
p
n
+
1
(
x
)
:=
a
n
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
n
)
{\displaystyle p_{n+1}(x):=a_{n}(x-x_{0})(x-x_{1})\cdots (x-x_{n})}
mit beliebigem
a
n
∈
R
∖
{
0
}
{\displaystyle a_{n}\in \mathbb {R} \setminus \{0\}}
.
Für
p
∈
Π
n
{\displaystyle p\in \Pi _{n}}
folgt
p
p
n
+
1
∈
Π
2
n
+
1
{\displaystyle pp_{n+1}\in \Pi _{2n+1}}
und somit
⟨
p
,
p
n
+
1
⟩
=
∫
a
b
w
(
x
)
p
(
x
)
p
n
+
1
(
x
)
d
x
=
I
n
(
p
p
n
+
1
)
=
∑
i
=
0
n
σ
i
p
(
x
i
)
p
n
+
1
(
x
i
)
⏟
=
0
=
0.
{\displaystyle \langle p,p_{n+1}\rangle =\int \limits _{a}^{b}w(x)p(x)p_{n+1}(x)\,dx={\mathcal {I}}_{n}(pp_{n+1})=\sum _{i=0}^{n}\sigma _{i}p(x_{i})\underbrace {p_{n+1}(x_{i})} _{=0}=0.}
q.e.d.
Zwei Funktionen
f
{\displaystyle f}
und
g
{\displaystyle g}
, für die
⟨
f
,
g
⟩
=
0
{\displaystyle \langle f,g\rangle =0}
gilt, nennt man orthogonal zueinander. Für eine Quadraturformel mit Genauigkeitsgrad
2
n
+
1
{\displaystyle 2n+1}
sollten die Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
also gerade als Nullstellen eines Polynoms vom Grad
n
+
1
{\displaystyle n+1}
gewählt werden, welches bezüglich des Skalarproduktes
⟨
⋅
,
⋅
⟩
{\displaystyle \langle \cdot ,\cdot \rangle }
orthogonal zu dem ganzen Raum
Π
n
{\displaystyle \Pi _{n}}
ist. Offenbar kann man ein solches Polynom gewinnen, indem man durch Anwendung des Gram-Schmidt-Orthogonalisierungsverfahren auf die Monome
1
,
x
,
…
,
x
n
+
1
{\displaystyle 1,x,\ldots ,x^{n+1}}
orthogonale Polynome
p
j
∈
Π
j
{\displaystyle p_{j}\in \Pi _{j}}
hinsichtlich
⟨
⋅
,
⋅
⟩
{\displaystyle \langle \cdot ,\cdot \rangle }
erzeugt. Diese Polynome haben nämlich die Eigenschaft
Π
k
=
span
{
p
0
,
p
1
,
…
,
p
k
}
,
⟨
p
i
,
p
j
⟩
=
δ
i
j
‖
p
i
‖
2
(
i
,
j
,
k
=
0
,
1
,
…
,
n
+
1
)
,
{\displaystyle \Pi _{k}=\operatorname {span} \{p_{0},p_{1},\ldots ,p_{k}\},\quad \langle p_{i},p_{j}\rangle =\delta _{ij}\|p_{i}\|^{2}\quad (i,j,k=0,1,\ldots ,n+1),}
so dass sich insbesondere jedes
p
∈
Π
n
{\displaystyle p\in \Pi _{n}}
mit gewissen
a
j
{\displaystyle a_{j}}
(
j
=
0
,
1
,
…
,
n
)
{\displaystyle (j=0,1,\ldots ,n)}
in der Form
p
(
x
)
=
∑
j
=
0
n
a
j
p
j
(
x
)
{\displaystyle p(x)=\sum _{j=0}^{n}a_{j}p_{j}(x)}
schreiben lässt und folglich mit den zugehörigen
a
j
{\displaystyle a_{j}}
gilt:
(8.42)
⟨
p
,
p
n
+
1
⟩
=
⟨
∑
j
=
0
n
a
j
p
j
(
x
)
,
p
n
+
1
⟩
=
∑
j
=
0
n
a
j
⟨
p
j
,
p
n
+
1
⟩
=
0
,
p
∈
Π
n
.
{\displaystyle \langle p,p_{n+1}\rangle =\left\langle \sum _{j=0}^{n}a_{j}p_{j}(x),p_{n+1}\right\rangle =\sum _{j=0}^{n}a_{j}\langle p_{j},p_{n+1}\rangle =0,\quad p\in \Pi _{n}.}
Darüber hinaus haben diese orthogonalen Polynome
p
j
{\displaystyle p_{j}}
nur reelle und einfache Nullstellen, welche alle im Intervall
(
a
,
b
)
{\displaystyle (a,b)}
liegen, wie im nächsten Unterabschnitt gezeigt wird. Die Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
sollten demzufolge gerade die Nullstellen des
(
n
+
1
)
{\displaystyle (n+1)}
-ten dieser orthogonalen Polynome sein. Die Gewichte
σ
i
{\displaystyle \sigma _{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
einer derartigen Quadraturformel sind dann gemäß Satz 8.4, der genauso für gewichtete Integrale gilt, durch
σ
i
:=
I
(
L
i
)
=
∫
a
b
w
(
x
)
L
i
(
x
)
d
x
=
⟨
L
i
,
1
⟩
,
i
=
0
,
1
,
…
,
n
{\displaystyle \sigma _{i}:=I(L_{i})=\int \limits _{a}^{b}w(x)L_{i}(x)\,dx=\langle L_{i},\mathbf {1} \rangle ,\quad i=0,1,\ldots ,n}
festgelegt, wobei
L
i
{\displaystyle L_{i}}
wieder die Lagrangeschen Basispolynome zu den Stützstellen
x
k
{\displaystyle x_{k}}
sind. Nach Definition 8.8 (entsprechend für gewichtete Integrale formuliert) handelt es sich bei der so definierten Formel um eine interpolatorische Quadraturformel.
Bevor wir diese sog. Gaußschen Quadraturformeln noch etwas näher betrachten, wollen wir auf ihren wesentlichen Baustein, orthogonale Polynome, näher eingehen.
Wie bereits im vorigen Unterabschnitt gesagt wurde, erhält man eine spezielle Folge
p
~
n
{\displaystyle {\tilde {p}}_{n}}
paarweise orthonormaler Polynome
p
~
n
∈
Π
n
{\displaystyle {\tilde {p}}_{n}\in \Pi _{n}}
durch Gram-Schmidt-Orthonormalisierung der Monome
1
,
x
,
x
2
,
…
{\displaystyle 1,x,x^{2},\ldots }
:
p
0
:=
1
,
p
~
0
:=
p
0
‖
p
0
‖
,
p
n
:=
x
n
−
∑
j
=
0
n
−
1
⟨
x
n
,
p
~
j
⟩
p
~
j
,
p
~
n
:=
p
n
‖
p
n
‖
(
n
=
1
,
2
,
…
)
.
{\displaystyle p_{0}:=1,\quad {\tilde {p}}_{0}:={\frac {p_{0}}{\|p_{0}\|}},\quad p_{n}:=x^{n}-\sum _{j=0}^{n-1}\langle x^{n},{\tilde {p}}_{j}\rangle {\tilde {p}}_{j},\quad {\tilde {p}}_{n}:={\frac {p_{n}}{\|p_{n}\|}}\quad (n=1,2,\ldots ).}
Statt mit den orthonormalen Polynomen
p
~
n
{\displaystyle {\tilde {p}}_{n}}
zu arbeiten, deren Hauptkoeffizienten i. a. von 1 verschieden sind, ist es manchmal bequemer, dies mit den orthogonalen Polynomen
p
n
∈
Π
n
{\displaystyle p_{n}\in \Pi _{n}}
zu tun, d. h. mit
(8.43)
p
0
:=
1
,
p
n
:=
x
n
−
∑
j
=
0
n
−
1
⟨
x
n
,
p
j
⟩
‖
p
j
‖
2
p
j
(
n
=
1
,
2
,
…
)
.
{\displaystyle p_{0}:=1,\quad p_{n}:=x^{n}-\sum _{j=0}^{n-1}{\frac {\langle x^{n},p_{j}\rangle }{\|p_{j}\|^{2}}}p_{j}\quad (n=1,2,\ldots ).}
Diese unterscheiden sich von den
p
~
n
{\displaystyle {\tilde {p}}_{n}}
nur durch den Skalar
1
/
‖
p
n
‖
{\displaystyle 1/\|p_{n}\|}
und haben offensichtlich den Hauptkoeffizienten 1. Für sie gilt
(8.44)
Π
k
=
span
{
p
0
,
p
1
,
…
,
p
k
}
,
⟨
p
i
,
p
j
⟩
=
δ
i
j
‖
p
i
‖
2
(
i
,
j
,
k
∈
N
0
)
{\displaystyle \Pi _{k}=\operatorname {span} \{p_{0},p_{1},\ldots ,p_{k}\},\quad \langle p_{i},p_{j}\rangle =\delta _{ij}\|p_{i}\|^{2}\quad (i,j,k\in \mathbb {N} _{0})}
und somit (vgl. (8.42))
(8.45)
⟨
p
,
p
k
⟩
=
0
,
p
∈
Π
k
−
1
(
k
=
1
,
2
,
…
)
.
{\displaystyle \langle p,p_{k}\rangle =0,\quad p\in \Pi _{k-1}\quad (k=1,2,\ldots ).}
Nach Konstruktion ist also
p
n
{\displaystyle p_{n}}
ein Polynom vom genauen Grad
n
{\displaystyle n}
mit Hauptkoeffizienten 1.
Die Polynome
p
n
{\displaystyle p_{n}}
können statt über die Formel (8.43) auch nach der im folgenden Satz angegebenen Drei-Term-Rekursionsformel berechnet werden.
Die Orthogonalpolynome in (8.43) genügen der Drei-Term-Rekursionsformel
p
0
(
x
)
=
1
,
p
1
(
x
)
=
x
−
β
0
,
p
n
+
1
(
x
)
=
(
x
−
β
n
)
p
n
(
x
)
−
γ
n
2
p
n
−
1
(
x
)
,
n
=
1
,
2
,
…
{\displaystyle p_{0}(x)=1,\quad p_{1}(x)=x-\beta _{0},\quad p_{n+1}(x)=(x-\beta _{n})p_{n}(x)-\gamma _{n}^{2}p_{n-1}(x),\quad n=1,2,\ldots }
mit den Koeffizienten
β
n
:=
⟨
x
p
n
,
p
n
⟩
‖
p
n
‖
2
(
n
=
0
,
1
,
…
)
,
γ
n
2
:=
‖
p
n
‖
2
‖
p
n
−
1
‖
2
(
n
=
1
,
2
,
…
)
.
{\displaystyle \beta _{n}:={\frac {\langle xp_{n},p_{n}\rangle }{\|p_{n}\|^{2}}}\quad (n=0,1,\ldots ),\quad \gamma _{n}^{2}:={\frac {\|p_{n}\|^{2}}{\|p_{n-1}\|^{2}}}\quad (n=1,2,\ldots ).}
Offenbar ist die behauptete Darstellung richtig für
p
0
{\displaystyle p_{0}}
und
p
1
{\displaystyle p_{1}}
. Für
n
≥
1
{\displaystyle n\geq 1}
setzen wir
q
n
+
1
:=
(
x
−
β
n
)
p
n
−
γ
n
2
p
n
−
1
{\displaystyle q_{n+1}:=(x-\beta _{n})p_{n}-\gamma _{n}^{2}p_{n-1}}
und zeigen im Folgenden
q
n
+
1
=
p
n
+
1
{\displaystyle q_{n+1}=p_{n+1}}
. Dazu stellen wir fest, dass
p
n
+
1
{\displaystyle p_{n+1}}
sowie
q
n
+
1
{\displaystyle q_{n+1}}
Polynome vom genauen Grad
n
+
1
{\displaystyle n+1}
sind und beide den Hauptkoeffizienten 1 haben. Somit gilt
(8.46)
r
:=
p
n
+
1
−
q
n
+
1
∈
Π
n
.
{\displaystyle r:=p_{n+1}-q_{n+1}\in \Pi _{n}.}
Wir zeigen nun, dass
q
n
+
1
{\displaystyle q_{n+1}}
wie
p
n
+
1
{\displaystyle p_{n+1}}
orthogonal zu dem ganzen Raum
Π
n
{\displaystyle \Pi _{n}}
ist und damit
(8.47)
⟨
p
,
q
n
+
1
⟩
=
0
,
p
∈
Π
n
{\displaystyle \langle p,q_{n+1}\rangle =0,\quad p\in \Pi _{n}}
gilt. Die Beziehungen (8.46) und (8.47) zusammen ergeben dann
‖
r
‖
2
=
⟨
r
,
r
⟩
=
⟨
r
,
p
n
+
1
−
q
n
+
1
⟨
=
⟨
r
,
p
n
+
1
⟩
−
⟨
r
,
q
n
+
1
⟩
=
0
{\displaystyle \|r\|^{2}=\langle r,r\rangle =\langle r,p_{n+1}-q_{n+1}\langle =\langle r,p_{n+1}\rangle -\langle r,q_{n+1}\rangle =0}
und folglich
r
=
0
{\displaystyle r=0}
bzw., wie behauptet,
p
n
+
1
=
q
n
+
1
{\displaystyle p_{n+1}=q_{n+1}}
.
Wir wollen nun (8.47) nachweisen. Aufgrund von
⟨
p
n
,
p
n
−
1
⟩
=
0
{\displaystyle \langle p_{n},p_{n-1}\rangle =0}
erhalten wir mit der Definition von
β
n
{\displaystyle \beta _{n}}
(8.48)
⟨
q
n
+
1
,
p
n
⟩
=
⟨
(
x
−
β
n
)
p
n
−
γ
n
2
p
n
−
1
,
p
n
⟩
=
⟨
x
p
n
,
p
n
⟩
−
β
n
‖
p
n
‖
2
=
0
{\displaystyle \langle q_{n+1},p_{n}\rangle =\langle (x-\beta _{n})p_{n}-\gamma _{n}^{2}p_{n-1},p_{n}\rangle =\langle xp_{n},p_{n}\rangle -\beta _{n}\|p_{n}\|^{2}=0}
und mit der Definition von
γ
n
{\displaystyle \gamma _{n}}
⟨
q
n
+
1
,
p
n
−
1
⟩
=
⟨
(
x
−
β
n
)
p
n
−
γ
n
2
p
n
−
1
,
p
n
−
1
⟩
=
⟨
x
p
n
,
p
n
−
1
⟩
−
γ
n
2
‖
p
n
−
1
‖
2
=
⟨
p
n
,
x
p
n
−
1
⟩
−
⟨
p
n
,
p
n
⟩
{\displaystyle \langle q_{n+1},p_{n-1}\rangle =\langle (x-\beta _{n})p_{n}-\gamma _{n}^{2}p_{n-1},p_{n-1}\rangle =\langle xp_{n},p_{n-1}\rangle -\gamma _{n}^{2}\|p_{n-1}\|^{2}=\langle p_{n},xp_{n-1}\rangle -\langle p_{n},p_{n}\rangle }
(8.49)
=
⟨
p
n
,
x
p
n
−
1
−
p
n
⟩
=
0
,
{\displaystyle =\langle p_{n},xp_{n-1}-p_{n}\rangle =0,}
wobei das letzte Gleichheitszeichen wegen
x
p
n
−
1
−
p
n
∈
Π
n
−
1
{\displaystyle xp_{n-1}-p_{n}\in \Pi _{n-1}}
gilt. Schließlich folgt:
(8.50)
⟨
q
n
+
1
,
p
j
⟩
=
⟨
p
n
,
x
p
j
⟩
⏟
=
0
−
β
n
⟨
p
n
,
p
j
⟩
⏟
=
0
−
γ
n
2
⟨
p
n
−
1
,
p
j
⟩
⏟
=
0
=
0
,
j
=
0
,
1
,
…
,
n
−
2.
{\displaystyle \langle q_{n+1},p_{j}\rangle =\underbrace {\langle p_{n},xp_{j}\rangle } _{=0}-\beta _{n}\underbrace {\langle p_{n},p_{j}\rangle } _{=0}-\gamma _{n}^{2}\underbrace {\langle p_{n-1},p_{j}\rangle } _{=0}=0,\quad j=0,1,\ldots ,n-2.}
Da sich jedes
p
∈
Π
n
{\displaystyle p\in \Pi _{n}}
gemäß (8.44) als Linearkombination der
p
j
{\displaystyle p_{j}}
(
j
=
0
,
1
,
…
,
n
)
{\displaystyle (j=0,1,\ldots ,n)}
darstellen lässt, folgt aus (8.48), (8.49) und (8.50) für jedes
p
∈
Π
n
{\displaystyle p\in \Pi _{n}}
mit gewissen
a
j
{\displaystyle a_{j}}
⟨
q
n
+
1
,
p
⟩
=
⟨
q
n
+
1
,
∑
j
=
0
n
a
j
p
j
(
x
)
⟩
=
∑
j
=
0
n
a
j
⟨
q
n
+
1
,
p
j
⟩
=
0.
{\displaystyle \langle q_{n+1},p\rangle =\left\langle q_{n+1},\sum _{j=0}^{n}a_{j}p_{j}(x)\right\rangle =\sum _{j=0}^{n}a_{j}\langle q_{n+1},p_{j}\rangle =0.}
Damit ist alles gezeigt.
q.e.d.
Für die Nullstellen der
p
j
{\displaystyle p_{j}}
(
j
∈
N
)
{\displaystyle (j\in \mathbb {N} )}
in (8.43) hat man folgende Aussage:
Die Nullstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
−
1
)
{\displaystyle (i=0,1,\ldots ,n-1)}
des
n
{\displaystyle n}
-ten Orthogonalpolynoms
p
n
{\displaystyle p_{n}}
in (8.43) sind reell, einfach und liegen alle in
(
a
,
b
)
{\displaystyle (a,b)}
. Sie besitzen die Darstellung
(8.51)
x
i
=
⟨
x
L
i
,
L
i
⟩
‖
L
i
‖
2
(
i
=
0
,
1
,
…
,
n
−
1
)
,
{\displaystyle x_{i}={\frac {\langle xL_{i},L_{i}\rangle }{\|L_{i}\|^{2}}}\quad (i=0,1,\ldots ,n-1),}
wobei
L
i
∈
Π
n
−
1
{\displaystyle L_{i}\in \Pi _{n-1}}
die zu den
x
k
{\displaystyle x_{k}}
(
k
=
0
,
1
,
…
,
n
−
1
)
{\displaystyle (k=0,1,\ldots ,n-1)}
gehörenden Lagrangeschen Basispolynome sind.
Es seien die Nullstellen
x
i
{\displaystyle x_{i}}
von
p
n
{\displaystyle p_{n}}
so durchnumeriert, dass
a
<
x
0
<
…
<
x
j
−
1
<
b
{\displaystyle a<x_{0}<\ldots <x_{j-1}<b}
(
0
≤
j
≤
n
)
{\displaystyle (0\leq j\leq n)}
diejenigen Nullstellen von
p
n
{\displaystyle p_{n}}
in
(
a
,
b
)
{\displaystyle (a,b)}
seien, an denen
p
n
{\displaystyle p_{n}}
sein Vorzeichen wechselt und die somit eine ungerade Vielfachheit haben. Wäre nun
j
−
1
≤
n
−
2
{\displaystyle j-1\leq n-2}
bzw.
j
≤
n
−
1
{\displaystyle j\leq n-1}
, so hätte das Polynom
q
(
x
)
:=
∏
k
=
0
j
−
1
(
x
−
x
k
)
{\displaystyle q(x):=\prod _{k=0}^{j-1}(x-x_{k})}
den Grad
j
≤
n
−
1
{\displaystyle j\leq n-1}
, so dass wegen (8.45)
(8.52)
⟨
p
n
,
q
⟩
=
0
{\displaystyle \langle p_{n},q\rangle =0}
folgte. Da die
x
k
{\displaystyle x_{k}}
(
k
=
j
,
j
+
1
,
…
,
n
−
1
)
{\displaystyle (k=j,j+1,\ldots ,n-1)}
Nullstellen von
p
n
{\displaystyle p_{n}}
mit gerader Vielfachheit wären, wäre dann aber
p
n
(
x
)
q
(
x
)
=
(
∏
k
=
0
j
−
1
(
x
−
x
k
)
2
)
(
∏
k
=
j
n
−
1
(
x
−
x
k
)
)
≥
0
,
x
∈
[
a
,
b
]
{\displaystyle p_{n}(x)q(x)=\left(\prod _{k=0}^{j-1}(x-x_{k})^{2}\right)\left(\prod _{k=j}^{n-1}(x-x_{k})\right)\geq 0,\quad x\in [a,b]}
und demzufolge
⟨
p
n
,
q
⟩
=
∫
a
b
w
(
x
)
p
n
(
x
)
q
(
x
)
d
x
>
0
{\displaystyle \langle p_{n},q\rangle =\int \limits _{a}^{b}w(x)p_{n}(x)q(x)\,dx>0}
im Widerspruch zu (8.52). Also ist
j
=
n
{\displaystyle j=n}
.
Um zur Darstellung (8.51) zu gelangen, schreibt man
p
n
{\displaystyle p_{n}}
für
n
≥
1
{\displaystyle n\geq 1}
in der Form
p
n
(
x
)
=
(
x
−
x
i
)
q
^
(
x
)
,
q
^
(
x
)
:=
∏
k
=
0
k
≠
i
n
−
1
(
x
−
x
k
)
.
{\displaystyle p_{n}(x)=(x-x_{i}){\hat {q}}(x),\quad {\hat {q}}(x):=\prod _{k=0 \atop k\neq i}^{n-1}(x-x_{k}).}
Es folgt
q
^
∈
Π
n
−
1
{\displaystyle {\hat {q}}\in \Pi _{n-1}}
sowie
0
=
⟨
p
n
,
q
^
⟩
=
⟨
x
q
^
,
q
^
⟨
−
x
i
⟨
q
^
,
q
^
⟩
.
{\displaystyle 0=\langle p_{n},{\hat {q}}\rangle =\langle x{\hat {q}},{\hat {q}}\langle -x_{i}\langle {\hat {q}},{\hat {q}}\rangle .}
Daraus ergibt sich wegen
⟨
q
^
,
q
^
⟩
≠
0
{\displaystyle \langle {\hat {q}},{\hat {q}}\rangle \neq 0}
x
i
=
⟨
x
q
^
,
q
^
⟩
‖
q
^
‖
2
=
⟨
x
L
i
,
L
i
⟩
‖
L
i
‖
2
{\displaystyle x_{i}={\frac {\langle x{\hat {q}},{\hat {q}}\rangle }{\|{\hat {q}}\|^{2}}}={\frac {\langle xL_{i},L_{i}\rangle }{\|L_{i}\|^{2}}}}
wobei sich die letzte Gleichung aus der Tatsache ergibt, dass das Polynom
q
^
{\displaystyle {\hat {q}}}
bis auf einen konstanten Faktor mit
L
i
{\displaystyle L_{i}}
übereinstimmt.
q.e.d.
In folgender Tabelle sind für verschiedene Intervalle und Gewichtsfunktionen die Bezeichnungen der zugehörigen orthogonalen Polynome aufgelistet:
Fehler beim Parsen (Konvertierungsfehler. Der Server („https://wikimedia.org/api/rest_“) hat berichtet: „Cannot get mml. TeX parse error: Bracket argument to \\ must be a dimension“): {\displaystyle {\begin{array}{|c|c|c|}\hline {\text{Intervall}}&{\text{Gewichtsfunktion }}w(x)&{\text{Name}}\\\hline [-1,1]&1&{\text{Legendre-Polynome}}\\[-1,1]&(1-x)^{\alpha }(1+x)^{\beta },\alpha ,\beta >-1&{\text{Jacobi-Polynome}}\\[-1,1]&1/{\sqrt {1-x^{2}}}&{\text{Tschebyscheff-Pol. der 1. Art}}\\[-1,1]&{\sqrt {1-x^{2}}}&{\text{Tschebyscheff-Pol. der 2. Art}}\\[0,\infty )&e^{-x^{2}}x^{\alpha },\alpha >-1&{\text{Laguerre-Polynome}}\\(-\infty ,\infty )&e^{-x^{2}}&{\text{Hermite-Polynome}}\\\hline \end{array}}}
Man kann zeigen (siehe z. B. E. W. Cheney: Introduction to Approximation Theory, 2nd ed., Chelsea Publish. Comp., New York, 1982):
Für
p
n
{\displaystyle p_{n}}
aus (8.43) gilt
‖
p
n
‖
=
min
a
0
,
…
,
a
n
−
1
∈
R
‖
x
n
+
a
n
−
1
x
n
−
1
+
…
+
a
1
x
+
a
0
‖
.
{\displaystyle \|p_{n}\|=\min _{a_{0},\ldots ,a_{n-1}\in \mathbb {R} }\left\|x^{n}+a_{n-1}x^{n-1}+\ldots +a_{1}x+a_{0}\right\|.}
Unter allen Polynomen vom Grad
n
{\displaystyle n}
mit Hauptkoeffizientem 1 macht also
p
n
{\displaystyle p_{n}}
die Norm in (8.40) minimal. Im Fall der Tschebyscheff-Polynome 1. Art minimiert
p
n
{\displaystyle p_{n}}
unter all diesen Polynomen überdies die Maximum-Norm (Satz 6.19) und im Fall der Tschebyscheff-Polynome 2. Art (s. Cheney) die (ungewichtete)
L
1
{\displaystyle L_{1}}
-Norm
‖
f
‖
:=
∫
a
b
|
f
(
x
)
|
d
x
,
f
∈
Π
.
{\displaystyle \|f\|:=\int \limits _{a}^{b}|f(x)|\,dx,\quad f\in \Pi .}
Mit Satz 8.30 sollen die Legendre-Polynome für
n
=
0
,
1
,
2
,
3
{\displaystyle n=0,1,2,3}
berechnet werden. Es ist somit
a
:=
−
1
,
b
:=
1
,
w
≡
1
{\displaystyle a:=-1,b:=1,w\equiv 1}
und folglich
β
n
=
⟨
x
p
n
,
p
n
⟩
‖
p
n
‖
2
=
∫
−
1
1
x
p
n
2
(
x
)
d
x
∫
−
1
1
p
n
2
(
x
)
d
x
,
γ
n
2
=
‖
p
n
‖
2
‖
p
n
−
1
‖
2
=
∫
−
1
1
p
n
2
(
x
)
d
x
∫
−
1
1
p
n
−
1
2
(
x
)
d
x
.
{\displaystyle \beta _{n}={\frac {\langle xp_{n},p_{n}\rangle }{\|p_{n}\|^{2}}}={\frac {\int _{-1}^{1}xp_{n}^{2}(x)\,dx}{\int _{-1}^{1}p_{n}^{2}(x)\,dx}},\quad \gamma _{n}^{2}={\frac {\|p_{n}\|^{2}}{\|p_{n-1}\|^{2}}}={\frac {\int _{-1}^{1}p_{n}^{2}(x)\,dx}{\int _{-1}^{1}p_{n-1}^{2}(x)\,dx}}.}
Mit
p
0
(
x
)
=
1
{\displaystyle p_{0}(x)=1}
ist
β
0
=
∫
−
1
1
x
d
x
∫
−
1
1
d
x
=
0
{\displaystyle \beta _{0}={\frac {\int _{-1}^{1}x\,dx}{\int _{-1}^{1}dx}}=0}
und damit weiter
p
1
(
x
)
=
x
{\displaystyle p_{1}(x)=x}
. Es ergeben sich ferner
β
1
=
∫
−
1
1
x
⋅
x
2
d
x
∫
−
1
1
x
2
d
x
=
0
,
γ
1
2
=
∫
−
1
1
x
2
d
x
∫
−
1
1
d
x
=
1
3
{\displaystyle \beta _{1}={\frac {\int _{-1}^{1}x\cdot x^{2}\,dx}{\int _{-1}^{1}x^{2}\,dx}}=0,\quad \gamma _{1}^{2}={\frac {\int _{-1}^{1}x^{2}\,dx}{\int _{-1}^{1}dx}}={\frac {1}{3}}}
und demnach
p
2
(
x
)
=
x
2
−
1
3
{\displaystyle p_{2}(x)=x^{2}-{\frac {1}{3}}}
sowie
β
2
=
∫
−
1
1
x
(
x
2
−
1
3
)
d
x
∫
−
1
1
(
x
2
−
1
3
)
2
d
x
=
0
,
γ
1
2
=
∫
−
1
1
(
x
2
−
1
3
)
2
d
x
∫
−
1
1
x
2
d
x
=
4
15
,
{\displaystyle \beta _{2}={\frac {\int _{-1}^{1}x\left(x^{2}-{\frac {1}{3}}\right)\,dx}{\int _{-1}^{1}\left(x^{2}-{\frac {1}{3}}\right)^{2}\,dx}}=0,\quad \gamma _{1}^{2}={\frac {\int _{-1}^{1}\left(x^{2}-{\frac {1}{3}}\right)^{2}\,dx}{\int _{-1}^{1}x^{2}\,dx}}={\frac {4}{15}},}
so dass folgt
p
3
(
x
)
=
x
3
−
1
3
x
−
4
15
x
=
x
3
−
3
5
x
.
{\displaystyle p_{3}(x)=x^{3}-{\frac {1}{3}}x-{\frac {4}{15}}x=x^{3}-{\frac {3}{5}}x.}
Für
n
∈
N
{\displaystyle n\in \mathbb {N} }
seien
p
j
{\displaystyle p_{j}}
(
j
=
0
,
1
,
…
,
n
,
n
+
1
)
{\displaystyle (j=0,1,\ldots ,n,n+1)}
die durch (8.43) definierten bezüglich
⟨
⋅
,
⋅
⟩
{\displaystyle \langle \cdot ,\cdot \rangle }
orthogonalen Polynome,
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
die Nullstellen von
p
n
+
1
{\displaystyle p_{n+1}}
und
σ
i
{\displaystyle \sigma _{i}}
die durch
σ
i
:=
⟨
L
i
,
1
⟩
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle \sigma _{i}:=\langle L_{i},1\rangle \quad (i=0,1,\ldots ,n)}
definierten Gewichte. Dann ist die Quadraturformel
(8.53)
I
n
(
f
)
:=
∑
i
=
0
n
σ
i
f
(
x
i
)
{\displaystyle {\mathcal {I}}_{n}(f):=\sum _{i=0}^{n}\sigma _{i}f(x_{i})}
interpolatorisch und hat (exakt) den Genauigkeitsgrad
2
n
+
1
{\displaystyle 2n+1}
.
Nach Definition 8.8 (entsprechend für gewichtete Integrale formuliert) ist
I
n
{\displaystyle {\mathcal {I}}_{n}}
aufgrund der Wahl der Gewichte eine interpolatorische Quadraturformel. Nach Korollar 8.9 hat eine solche mindestens den Genauigkeitsgrad
n
{\displaystyle n}
. Wir wollen nun zeigen, dass sie mindestens den Genauigkeitsgrad
2
n
+
1
{\displaystyle 2n+1}
und damit exakt den Genauigkeitsgrad
2
n
+
1
{\displaystyle 2n+1}
besitzt, wie aus den Argumenten in Abschnitt 8.5.1 hervorgeht.
Es sei
p
∈
Π
2
n
+
1
{\displaystyle p\in \Pi _{2n+1}}
beliebig. Dann lässt sich
p
{\displaystyle p}
mit gewissen
q
,
r
∈
Π
n
{\displaystyle q,r\in \Pi _{n}}
nach einer Polynomdivision mit Rest in der Form
p
=
q
p
n
+
1
+
r
{\displaystyle p=qp_{n+1}+r}
schreiben. Wegen
p
n
+
1
(
x
i
)
=
0
{\displaystyle p_{n+1}(x_{i})=0}
gilt dann
p
(
x
i
)
=
r
(
x
i
)
,
i
=
0
,
1
,
…
,
n
.
{\displaystyle p(x_{i})=r(x_{i}),\quad i=0,1,\ldots ,n.}
Mit der Lagrangeschen Interpolationsformel (6.7), angewandt auf
r
{\displaystyle r}
, erhält man weiter
r
(
x
)
=
∑
i
=
0
n
r
(
x
i
)
L
i
(
x
)
=
∑
i
=
0
n
p
(
x
i
)
L
i
(
x
)
.
{\displaystyle r(x)=\sum _{i=0}^{n}r(x_{i})L_{i}(x)=\sum _{i=0}^{n}p(x_{i})L_{i}(x).}
Somit schließt man
(8.54)
I
(
p
)
=
⟨
p
,
1
⟩
=
⟨
q
,
p
n
+
1
⟩
⏟
=
0
+
⟨
r
,
1
⟩
=
∑
i
=
0
n
p
(
x
i
)
⟨
L
i
,
1
⟩
=
∑
i
=
0
n
σ
i
p
(
x
i
)
,
{\displaystyle {\mathcal {I}}(p)=\langle p,\mathbf {1} \rangle =\underbrace {\langle q,p_{n+1}\rangle } _{=0}+\langle r,\mathbf {1} \rangle =\sum _{i=0}^{n}p(x_{i})\langle L_{i},\mathbf {1} \rangle =\sum _{i=0}^{n}\sigma _{i}p(x_{i}),}
womit der Genauigkeitsgrad von mindestens
2
n
+
1
{\displaystyle 2n+1}
für
I
n
{\displaystyle {\mathcal {I}}_{n}}
nachgewiesen ist.
q.e.d.
Die Quadraturformel in (8.53) mit den in Satz 8.34 genannten Stützstellen
x
i
{\displaystyle x_{i}}
und Gewichten
σ
i
{\displaystyle \sigma _{i}}
bezeichnet man als Gaußsche Quadraturformel . Ihr Genauigkeitsgrad ist optimal, da es keine Quadraturformeln mit Genauigkeitsgrad
2
n
+
2
{\displaystyle 2n+2}
gibt (vgl. Abschnitt 8.5.1). Weiter hat man:
Für die Gewichte
σ
i
{\displaystyle \sigma _{i}}
der Gaußschen Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
von Satz 8.34 gilt
σ
i
=
⟨
L
i
,
L
i
⟩
>
0
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle \sigma _{i}=\langle L_{i},L_{i}\rangle >0\quad (i=0,1,\ldots ,n)}
und
(8.55)
∑
i
=
0
n
σ
i
=
∫
a
b
w
(
x
)
d
x
.
{\displaystyle \sum _{i=0}^{n}\sigma _{i}=\int \limits _{a}^{b}w(x)\,dx.}
Wendet man die Beziehungen (8.54) auf
p
:=
L
j
2
∈
Π
2
n
{\displaystyle p:=L_{j}^{2}\in \Pi _{2n}}
an, so folgt
⟨
L
j
,
L
j
⟩
=
⟨
L
j
2
,
1
⟩
=
∑
i
=
0
n
σ
i
L
j
2
(
x
i
)
=
σ
j
.
{\displaystyle \langle L_{j},L_{j}\rangle =\left\langle L_{j}^{2},\mathbf {1} \right\rangle =\sum _{i=0}^{n}\sigma _{i}L_{j}^{2}(x_{i})=\sigma _{j}.}
Weiter gilt
⟨
L
j
,
L
j
⟩
>
0
,
d
a
m
a
n
<
m
a
t
h
>
L
j
2
(
x
)
≥
0
,
x
∈
[
a
,
b
]
{\displaystyle \langle L_{j},L_{j}\rangle >0,daman<math>L_{j}^{2}(x)\geq 0,x\in [a,b]}
sowie
L
j
2
(
x
)
>
0
{\displaystyle L_{j}^{2}(x)>0}
z. B. für alle
x
∈
(
x
j
−
1
,
x
j
+
1
)
{\displaystyle x\in (x_{j-1},x_{j+1})}
hat. Die Beziehung (8.55) folgt nun mit Satz 8.34 aus
∑
i
=
0
n
σ
i
=
I
n
(
1
)
=
I
(
1
)
=
∫
a
b
w
(
x
)
d
x
.
{\displaystyle \sum _{i=0}^{n}\sigma _{i}={\mathcal {I}}_{n}(1)={\mathcal {I}}(1)=\int \limits _{a}^{b}w(x)\,dx.}
q.e.d.
Anders als bei den abgeschlossenen Newton-Cotes-Formeln haben also die Gaußschen Quadraturformeln für alle
n
∈
N
0
{\displaystyle n\in \mathbb {N} _{0}}
positive Gewichte. Mit dem folgenden Satz geben wir abschließend eine Darstellung für den bei der Gauß-Quadratur entstehenden Quadraturfehler an.
Es sei
f
∈
C
2
n
+
2
[
a
,
b
]
{\displaystyle f\in C^{2n+2}[a,b]}
und In die Gaußsche Quadraturformel mit Stützstellen
x
k
{\displaystyle x_{k}}
(
k
=
0
,
1
,
…
,
n
)
{\displaystyle (k=0,1,\ldots ,n)}
. Dann gilt für
γ
^
2
n
+
1
:=
∫
0
1
w
(
(
b
−
a
)
t
+
a
)
2
∏
k
=
0
n
(
t
−
t
k
)
2
d
t
{\displaystyle {\hat {\gamma }}_{2n+1}:=\int \limits _{0}^{1}w((b-a)t+a)^{2}\prod _{k=0}^{n}(t-t_{k})^{2}\,dt}
mit
t
k
:=
x
k
−
a
b
−
a
(
k
=
0
,
1
,
…
,
n
)
{\displaystyle t_{k}:={\frac {x_{k}-a}{b-a}}\quad (k=0,1,\ldots ,n)}
und für ein
ξ
∈
[
a
,
b
]
{\displaystyle \xi \in [a,b]}
:
I
(
f
)
−
I
n
(
f
)
=
γ
^
2
n
+
1
(
b
−
a
)
2
n
+
3
(
2
n
+
2
)
!
f
(
2
n
+
2
)
(
ξ
)
.
{\displaystyle {\mathcal {I}}(f)-{\mathcal {I}}_{n}(f)={\hat {\gamma }}_{2n+1}{\frac {(b-a)^{2n+3}}{(2n+2)!}}f^{(2n+2)}(\xi ).}
Den Satz 8.15 kann man auf den Fall gewichteter Integrale übertragen und dann aufgrund von Satz 8.34 auf die Gaußsche Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
mit
r
:=
2
n
+
1
{\displaystyle r:=2n+1}
anwenden. Man wählt dort zu den Stützpunkten
x
k
{\displaystyle x_{k}}
von
I
n
{\displaystyle {\mathcal {I}}_{n}}
die weiteren Stützpunkte
t
n
+
1
:=
t
0
,
…
,
t
2
n
+
1
:=
t
n
{\displaystyle t_{n+1}:=t_{0},\ldots ,t_{2n+1}:=t_{n}}
, so dass insbesondere
∏
k
=
0
2
n
+
1
(
t
−
t
k
)
=
∏
k
=
0
n
(
t
−
t
k
)
2
=
p
n
+
1
2
(
x
)
≥
0
,
t
∈
[
0
,
1
]
{\displaystyle \prod _{k=0}^{2n+1}(t-t_{k})=\prod _{k=0}^{n}(t-t_{k})^{2}=p_{n+1}^{2}(x)\geq 0,\quad t\in [0,1]}
gilt. Weiter folge man dann dem Beweis von Satz 8.15.
q.e.d.
Natürlich ist es auch möglich, summierte Gaußsche Quadraturformeln zu definieren und zu verwenden. Die Resultate in Abschnitt 8.3 lassen sich ganz kanonisch auf solche Formeln übertragen.
8.5.4 Berechnung der Stützstellen und Gewichte
Bearbeiten
Beispielsweise für die Tschebyscheff-Polynome 1. Art kann man die Nullstellen explizit angeben (vgl. Satz 6.18). Im allgemeinen steht man bei Verwendung der Gaußschen Quadraturformeln für größere Werte von
n
{\displaystyle n}
aber vor dem Problem, die
n
+
1
{\displaystyle n+1}
Nullstellen
x
i
{\displaystyle x_{i}}
des Polynoms
p
n
+
1
{\displaystyle p_{n+1}}
der bezüglich
⟨
⋅
,
⋅
⟩
{\displaystyle \langle \cdot ,\cdot \rangle }
orthogonalen Polynome
p
j
{\displaystyle p_{j}}
(
j
=
0
,
1
,
2
,
…
)
{\displaystyle (j=0,1,2,\ldots )}
und/oder die Gewichte
σ
i
:=
⟨
L
i
,
1
⟩
{\displaystyle \sigma _{i}:=\langle L_{i},\mathbf {1} \rangle }
zu bestimmen. Wir wollen hier abschließend einen Weg zu ihrer Berechnung aufzeigen. Dazu gehen wir davon aus, dass die Koeffizienten
β
j
{\displaystyle \beta _{j}}
und
γ
j
{\displaystyle \gamma _{j}}
in der Rekursionsformel (8.43) für die orthogonalen
p
j
{\displaystyle p_{j}}
(8.56)
p
0
:=
1
,
p
1
:=
x
−
β
0
,
{\displaystyle p_{0}:=1,\quad p_{1}:=x-\beta _{0},}
(8.57)
p
k
+
1
:=
(
x
−
β
k
)
p
k
−
γ
k
2
p
k
−
1
,
k
=
1
,
2
,
…
{\displaystyle p_{k+1}:=(x-\beta _{k})pk-\gamma _{k}^{2}p_{k-1},\quad k=1,2,\ldots }
bereits bekannt sind und somit die symmetrische Tridiagonalmatrix
(8.58)
J
:=
(
β
0
−
γ
1
0
…
0
−
γ
1
β
1
−
γ
2
⋱
⋮
0
−
γ
2
⋱
⋱
0
⋮
⋱
⋱
⋱
−
γ
n
0
…
0
−
γ
n
β
n
)
∈
R
(
n
+
1
)
×
(
n
+
1
)
{\displaystyle J:={\begin{pmatrix}\beta _{0}&-\gamma _{1}&0&\ldots &0\\-\gamma _{1}&\beta _{1}&-\gamma _{2}&\ddots &\vdots \\0&-\gamma _{2}&\ddots &\ddots &0\\\vdots &\ddots &\ddots &\ddots &-\gamma _{n}\\0&\ldots &0&-\gamma _{n}&\beta _{n}\end{pmatrix}}\in \mathbb {R} ^{(n+1)\times (n+1)}}
aufgestellt werden kann. Der folgende Satz besagt nun, dass die Stützstellen
x
i
{\displaystyle x_{i}}
der Gaußschen Quadraturformeln die Eigenwerte von
J
{\displaystyle J}
sind und sich deren Gewichte
σ
i
{\displaystyle \sigma _{i}}
aus zugehörigen Eigenvektoren bestimmen lassen.
Für die Stützstellen
x
i
{\displaystyle x_{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
und die Gewichte
σ
i
{\displaystyle \sigma _{i}}
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle (i=0,1,\ldots ,n)}
der Gaußschen Quadraturformel
I
n
{\displaystyle {\mathcal {I}}_{n}}
gilt mit
v
(
i
)
:=
(
τ
0
p
0
(
x
i
)
⏟
=
1
,
…
,
τ
n
p
n
(
x
i
)
)
T
{\displaystyle v^{(i)}:=(\underbrace {\tau _{0}p_{0}(x_{i})} _{=1},\ldots ,\tau _{n}p_{n}(x_{i}))^{T}}
für
τ
k
:=
{
1
f
a
l
l
s
k
=
0
,
(
−
1
)
k
/
(
γ
1
γ
2
⋯
γ
k
)
f
a
l
l
s
k
∈
{
1
,
…
,
n
}
,
{\displaystyle \tau _{k}:={\begin{cases}1&falls\ k=0,\\(-1)^{k}/(\gamma _{1}\gamma _{2}\cdots \gamma _{k})&falls\ k\in \{1,\ldots ,n\},\end{cases}}}
wobei
p
k
{\displaystyle p_{k}}
(
k
=
0
,
1
,
…
,
n
)
{\displaystyle (k=0,1,\ldots ,n)}
die bezüglich
⟨
⋅
,
⋅
⟩
{\displaystyle \langle \cdot ,\cdot \rangle }
orthogonalen Polynome seien:
(8.59)
J
v
(
i
)
=
x
i
v
(
i
)
(
i
=
0
,
1
,
…
,
n
)
{\displaystyle Jv^{(i)}=x_{i}v^{(i)}\quad (i=0,1,\ldots ,n)}
und
(8.60)
σ
i
=
⟨
1
,
1
⟩
(
v
(
i
)
)
T
v
(
i
)
=
⟨
1
,
1
⟩
∑
k
=
0
n
τ
k
2
p
k
2
(
x
i
)
(
i
=
0
,
1
,
…
,
n
)
.
{\displaystyle \sigma _{i}={\frac {\langle \mathbf {1} ,\mathbf {1} \rangle }{(v^{(i)})^{T}v^{(i)}}}={\frac {\langle \mathbf {1} ,\mathbf {1} \rangle }{\sum _{k=0}^{n}\tau _{k}^{2}p_{k}^{2}(x_{i})}}\quad (i=0,1,\ldots ,n).}
Aus den Definitionen von
J
,
v
(
i
)
,
p
1
{\displaystyle J,v^{(i)},p_{1}}
und
τ
1
{\displaystyle \tau _{1}}
ergibt sich
(
J
v
(
i
)
)
1
=
β
0
−
γ
1
τ
1
p
1
(
x
i
)
=
β
0
+
p
1
(
x
i
)
=
β
0
+
x
i
β
0
=
x
i
=
x
i
v
1
(
i
)
.
{\displaystyle \left(Jv^{(i)}\right)_{1}=\beta _{0}-\gamma _{1}\tau _{1}p_{1}(x_{i})=\beta _{0}+p_{1}(x_{i})=\beta _{0}+x_{i}\beta _{0}=x_{i}=x_{i}v_{1}^{(i)}.}
Definiert man
γ
n
+
1
=
τ
n
+
1
:=
0
{\displaystyle \gamma _{n+1}=\tau _{n+1}:=0}
und berücksichtigt man
p
n
+
1
(
x
i
)
=
0
{\displaystyle p_{n+1}(x_{i})=0}
, so erhält man aus den Rekursionsformeln (8.56) und (8.57) mit
x
=
x
i
{\displaystyle x=x_{i}}
weiter
(
J
v
(
i
)
)
k
+
1
=
−
γ
k
τ
k
−
1
p
k
−
1
(
x
i
)
+
β
k
τ
k
p
k
(
x
i
)
−
γ
k
+
1
τ
k
+
1
p
k
+
1
(
x
i
)
{\displaystyle \left(Jv^{(i)}\right)_{k+1}=-\gamma _{k}\tau _{k-1}p_{k-1}(x_{i})+\beta _{k}\tau _{k}p_{k}(x_{i})-\gamma _{k+1}\tau _{k+1}p_{k+1}(x_{i})}
=
−
γ
k
(
−
γ
k
τ
k
)
p
k
−
1
(
x
i
)
+
β
k
τ
k
p
k
(
x
i
)
−
γ
k
+
1
(
−
1
)
γ
k
+
1
p
k
+
1
(
x
i
)
{\displaystyle =-\gamma _{k}(-\gamma _{k}\tau _{k})p_{k-1}(x_{i})+\beta _{k}\tau _{k}p_{k}(x_{i})-\gamma _{k+1}{\frac {(-1)}{\gamma _{k+1}}}p_{k+1}(x_{i})}
=
τ
k
[
γ
k
2
p
k
−
1
(
x
i
)
+
β
k
p
k
(
x
i
)
+
(
x
i
−
β
k
)
p
k
(
x
i
)
−
γ
k
2
p
k
−
1
(
x
i
)
]
=
τ
k
x
i
p
k
(
x
i
)
=
x
i
v
k
+
1
(
i
)
,
k
=
1
,
…
,
n
.
{\displaystyle =\tau _{k}\left[\gamma _{k}^{2}p_{k-1}(x_{i})+\beta _{k}p_{k}(x_{i})+(x_{i}-\beta _{k})p_{k}(x_{i})-\gamma _{k}^{2}p_{k-1}(x_{i})\right]=\tau _{k}x_{i}p_{k}(x_{i})=x_{i}v_{k+1}^{(i)},\quad k=1,\ldots ,n.}
Damit ist (8.59) bewiesen.
Gleichung (8.59) besagt, dass
v
(
i
)
{\displaystyle v^{(i)}}
Eigenvektor zum Eigenwert
x
i
{\displaystyle x_{i}}
der Matrix
J
{\displaystyle J}
ist. Gemäß Satz 8.31 sind diese Eigenwerte paarweise verschieden. Da für eine reelle symmetrische Matrix Eigenvektoren zu paarweise verschiedenen Eigenwerten orthogonal zueinander sind, folgt
(8.61)
(
v
(
i
)
)
T
v
(
k
)
=
0
(
i
≠
k
)
.
{\displaystyle \left(v^{(i)}\right)^{T}v^{(k)}=0\quad (i\neq k).}
Da ferner die Polynome
p
j
{\displaystyle p_{j}}
(
j
=
0
,
1
,
2
,
…
)
{\displaystyle (j=0,1,2,\ldots )}
paarweise orthogonal sind und die Gaußsche Quadraturformel alle
p
k
{\displaystyle p_{k}}
(
k
=
0
,
1
,
…
,
n
)
{\displaystyle (k=0,1,\ldots ,n)}
exakt integriert, folgt weiter mit Satz 8.34
(8.62)
δ
k
0
⟨
1
,
1
⟩
=
⟨
p
k
,
p
0
⟩
=
I
(
p
k
)
=
∑
j
=
0
n
σ
j
p
k
(
x
j
)
(
k
=
0
,
1
,
…
,
n
)
,
{\displaystyle \delta _{k0}\langle \mathbf {1} ,\mathbf {1} \rangle =\langle p_{k},p_{0}\rangle ={\mathcal {I}}(p_{k})=\sum _{j=0}^{n}\sigma _{j}p_{k}(x_{j})\quad (k=0,1,\ldots ,n),}
wobei
δ
i
j
{\displaystyle \delta _{ij}}
das Kroneckersymbol ist. Multiplikation von (8.62) mit
τ
k
2
p
k
(
x
i
)
{\displaystyle \tau _{k}^{2}p_{k}(x_{i})}
und anschließende Summation über
k
{\displaystyle k}
liefert schließlich unter Verwendung von (8.61)
⟨
1
,
1
⟩
∑
k
=
0
n
τ
k
2
p
k
(
x
i
)
δ
k
0
=
⟨
1
,
1
⟩
τ
0
2
p
0
(
x
i
)
=
⟨
1
,
1
⟩
=
∑
k
=
0
n
∑
j
=
0
n
σ
j
τ
k
2
p
k
(
x
i
)
p
k
(
x
j
)
{\displaystyle \langle \mathbf {1} ,\mathbf {1} \rangle \sum _{k=0}^{n}\tau _{k}^{2}p_{k}(x_{i})\delta _{k0}=\langle \mathbf {1} ,\mathbf {1} \rangle \tau _{0}^{2}p_{0}(x_{i})=\langle \mathbf {1} ,\mathbf {1} \rangle =\sum _{k=0}^{n}\sum _{j=0}^{n}\sigma _{j}\tau _{k}^{2}p_{k}(x_{i})p_{k}(x_{j})}
=
∑
j
=
0
n
σ
j
∑
k
=
0
n
τ
k
2
p
k
(
x
i
)
p
k
(
x
j
)
=
∑
j
=
0
n
σ
j
(
v
(
i
)
)
T
v
(
j
)
=
σ
i
(
v
(
i
)
)
T
v
(
i
)
{\displaystyle =\sum _{j=0}^{n}\sigma _{j}\sum _{k=0}^{n}\tau _{k}^{2}p_{k}(x_{i})p_{k}(x_{j})=\sum _{j=0}^{n}\sigma _{j}\left(v^{(i)}\right)^{T}v^{(j)}=\sigma _{i}\left(v^{(i)}\right)^{T}v^{(i)}}
Damit ist auch die Gültigkeit von (8.60) bewiesen.
q.e.d.
Wir verwenden Satz 8.37 für die Herleitung der Gaußschen Quadraturformel mit
n
:=
2
{\displaystyle n:=2}
zur näherungsweisen Berechnung des Integrals
∫
−
1
1
f
(
x
)
d
x
.
{\displaystyle \int \limits _{-1}^{1}f(x)\,dx.}
Beispiel 8.33 liefert
J
:=
(
0
−
1
/
3
0
−
1
/
3
0
−
2
/
15
0
−
2
/
15
0
)
.
{\displaystyle J:={\begin{pmatrix}0&-1/{\sqrt {3}}&0\\-1/{\sqrt {3}}&0&-2/{\sqrt {15}}\\0&-2/{\sqrt {15}}&0\end{pmatrix}}.}
Die Eigenwerte von
J
{\displaystyle J}
berechnen sich aus
det
(
J
−
λ
I
)
=
−
λ
(
λ
2
−
4
15
)
+
1
3
λ
=
λ
(
3
5
−
λ
2
)
{\displaystyle \det(J-\lambda I)=-\lambda \left(\lambda ^{2}-{\frac {4}{15}}\right)+{\frac {1}{3}}\lambda =\lambda \left({\frac {3}{5}}-\lambda ^{2}\right)}
so dass man die Stützstellen
x
0
=
−
3
/
5
,
x
1
=
0
{\displaystyle x_{0}=-{\sqrt {3/5}},x_{1}=0}
und
x
2
=
3
/
5
{\displaystyle x_{2}={\sqrt {3/5}}}
erhält. Weiter hat man
τ
0
=
1
,
τ
1
=
−
1
/
γ
1
=
−
3
,
τ
2
=
1
/
(
γ
1
γ
2
)
=
3
5
/
2
{\displaystyle \tau _{0}=1,\quad \tau _{1}=-1/\gamma _{1}=-{\sqrt {3}},\quad \tau _{2}=1/(\gamma _{1}\gamma _{2})=3{\sqrt {5}}/2}
sowie
⟨
1
,
1
⟩
=
∫
−
1
1
d
x
=
2.
{\displaystyle \langle \mathbf {1} ,\mathbf {1} \rangle =\int \limits _{-1}^{1}dx=2.}
Mit
p
0
(
x
)
=
1
,
p
1
(
x
)
=
x
,
p
2
(
x
)
=
x
2
−
1
3
{\displaystyle p_{0}(x)=1,\quad p_{1}(x)=x,\quad p_{2}(x)=x^{2}-{\frac {1}{3}}}
hat man
i
p
0
(
x
i
)
p
1
(
x
i
)
p
2
(
x
i
)
0
1
−
3
/
5
4
/
15
1
1
0
−
1
/
3
2
1
3
/
5
4
/
15
{\displaystyle {\begin{array}{|c||c|c|c|}\hline i&p_{0}(x_{i})&p_{1}(x_{i})&p_{2}(x_{i})\\\hline 0&1&-{\sqrt {3/5}}&4/15\\1&1&0&-1/3\\2&1&{\sqrt {3/5}}&4/15\\\hline \end{array}}}
und demnach
σ
0
=
2
τ
0
2
p
0
2
(
x
0
)
+
τ
1
2
p
1
2
(
x
0
)
+
τ
2
2
p
2
2
(
x
0
)
=
5
9
,
{\displaystyle \sigma _{0}={\frac {2}{\tau _{0}^{2}p_{0}^{2}(x_{0})+\tau _{1}^{2}p_{1}^{2}(x_{0})+\tau _{2}^{2}p_{2}^{2}(x_{0})}}={\frac {5}{9}},}
σ
1
=
2
τ
0
2
p
0
2
(
x
1
)
+
τ
1
2
p
1
2
(
x
1
)
+
τ
2
2
p
2
2
(
x
1
)
=
8
9
,
{\displaystyle \sigma _{1}={\frac {2}{\tau _{0}^{2}p_{0}^{2}(x_{1})+\tau _{1}^{2}p_{1}^{2}(x_{1})+\tau _{2}^{2}p_{2}^{2}(x_{1})}}={\frac {8}{9}},}
σ
2
=
2
τ
0
2
p
0
2
(
x
2
)
+
τ
1
2
p
1
2
(
x
2
)
+
τ
2
2
p
2
2
(
x
2
)
=
5
9
.
{\displaystyle \sigma _{2}={\frac {2}{\tau _{0}^{2}p_{0}^{2}(x_{2})+\tau _{1}^{2}p_{1}^{2}(x_{2})+\tau _{2}^{2}p_{2}^{2}(x_{2})}}={\frac {5}{9}}.}
Man erhält also die Gaußsche Quadraturformel
I
2
(
f
)
:=
5
9
f
(
−
3
5
)
+
8
9
f
(
0
)
+
5
9
f
(
3
5
)
{\displaystyle {\mathcal {I}}_{2}(f):={\frac {5}{9}}f\left(-{\sqrt {\frac {3}{5}}}\right)+{\frac {8}{9}}f(0)+{\frac {5}{9}}f\left({\sqrt {\frac {3}{5}}}\right)}
Die gesuchten Eigenwerte von
J
{\displaystyle J}
müssen aber, zumindest für größere
n
{\displaystyle n}
, normalerweise mit einer numerischen Methode bestimmt werden.