Kap4 Numeriche Mathe

Als pdf oder txt herunterladen
Als pdf oder txt herunterladen
Sie sind auf Seite 1von 50

Kapitel 4.

Numerische Integration

4.1 Interpolatorische Quadraturformeln


4.2 Gaußsche Quadraturformeln
4.3 Das Rombergsche Integrationsverfahren
4.4 Praktische Aspekte der Integration

Numerische Mathematik I 147


Interpolatorische Quadraturformeln

4.1 Interpolatorische Quadraturformeln

Ziel: Näherungswerte für bestimmte Integrale, wenn sich keine geschlossene Form
der Stammfunktion finden lässt.
Beispiele:
Z x
2
Die Stammfunktion F (x) = e −t dt benötigt man in der
0
Wahrscheinlichkeitstheorie, um Werte der Normalverteilung zu berechnen.
Z x
sin t
Der Integralsinus dt spielt eine Rolle in der Signalverarbeitung
0 t
Z x t Z x
e dt
Integral-Exponentialfunktion dt und Integral-Logarithmus
−∞ t 0 ln t

Numerische Mathematik I 148


Interpolatorische Quadraturformeln

Möglichkeiten zur numerischen Berechnung von Näherungswerten:


Potenzreihe der Stammfunktion bilden und die Partialsummen an den
Integrationsgrenzen auswerten

Hier: Formeln für Näherungsflächen verwenden (z.B. Trapezregel, Keplersche


Fass-Regel)

Numerische Mathematik I 149


Interpolatorische Quadraturformeln Definition: Quadraturformel

4.1.1 Definition: Quadraturformel


Z b
Eine Näherung des bestimmten Integrals f (x) dx wird durch die
a
Quadraturformel
n
X
In (f ; [a, b]) = ωi f (xi )
i =0

angegeben, wobei
x0 < x1 < · · · < xn die Knoten (meistens in [a, b] gewählt) und
ω0 , ω1 , . . . , ωn ∈ R die Gewichte
bezeichnen.
Die Quadraturformel heißt interpolatorisch, wenn ihre Gewichte
Z b
ωk = Ln,k (x) dx, k = 0, . . . , n
a

Y n
x − xj
sind, wobei Ln,k (x) = die Lagrange-Grundpolynome bezeichnen.
j =0
xk − xj
j 6=k

Numerische Mathematik I 150


Interpolatorische Quadraturformeln Definition: Quadraturformel

Beachte: Für jedes Polynom p ∈ Pn gilt


n
X
p(x) = p(xk )Ln,k (x). “p interpoliert sich selbst”
k=0

Also liefert jede interpolatorische Quadraturformel mit n + 1 Knoten


n
X Z b
In (p; [a, b]) = ωk p(xk ) = p(x) dx,
k=0 a

d.h. sie ergibt den exakten Wert des Integrals von p.

Definition und Bemerkung


Eine Quadraturformel hat den Exaktheitsgrad m (oder die Ordnung m + 1), falls
sie für alle Polynome p ∈ Pm vom Grad kleiner oder gleich m den exakten Wert
Rb
a p(x) dx liefert.

Jede interpolatorische Quadraturformel (mit n + 1 Knoten) hat mindestens den


Exaktheitsgrad n und höchstens den Exaktheitsgrad 2n + 1.

Denn: mindestens Exaktheitsgrad n: Definition interpolatorisch


Q R
höchstens Exaktheitsgrad 2n + 1: für q(x) = nj=0 (x − xj )2 gilt ab q(x) dx > 0, aber In (q) = 0.
Numerische Mathematik I 151
Interpolatorische Quadraturformeln Beispiele:

Z b
4.1.2 Beispiele: Für f (x) dx verwendete Quadraturformeln
a
 
a+b
a) (b − a)f Mittelpunktsregel
2
b−a
b) (f (a) + f (b)) Trapezregel
2
   
b−a a+b
c) f (a) + 4f + f (b) Simpsonregel
6 2
auch als Keplersche Fass-Regel bezeichnet.

Numerische Mathematik I 152


Interpolatorische Quadraturformeln Beispiele:

Durch Aufteilen des Intervalls [a, b] in N Teile [ak−1 , ak ] mit

b−a
ak = a + k · , k = 0, 1, . . . , N,
N
ak−1 +ak
und verwenden der Mittelpunkte yk = 2 , also

a = a0 < y1 < a1 < y2 < · · · < aN−1 < yN < aN = b

erhält man die “summierten” Regeln:


a’) summierte Mittelpunktsregel
b−a
I0Σ (f ; N) = (f (y1 ) + f (y2 ) + · · · + f (yN ))
N
b’) summierte Trapezregel
I1Σ (f ; N) = fracb − a2N (f (a) + 2f (a1 ) + 2f (a2 ) + · · · + 2f (aN−1 ) + f (b))

c’) summierte Simpsonregel


b−a
I2Σ (f ; N) = (f (a) + 4f (y1 ) + 2f (a1 ) + 4f (y2 ) + · · · + 2f (aN−1 ) + 4f (yN ) + f (b))
6N

Numerische Mathematik I 153


Interpolatorische Quadraturformeln Definition: Newton-Cotes Formeln

Systematik zur Definition der “einfachen” Regeln: Gegeben sei das Integral
Rb
a f (x) dx.

4.1.3 Definition: Newton-Cotes Formeln


Für gegebenes n ∈ N wählen wir
b−a
die Knoten xk = a + k · , k = 0, 1, . . . , n,
n
Z b
die Gewichte ωk = Ln,k (x) dx wie üblich für eine interpolatorische
a
Quadraturformel.

Bemerkung: Der Exaktheitsgrad der Newton-Cotes Formeln mit n + 1 Knoten ist


n für ungerades n, n + 1 für gerades n.

Numerische Mathematik I 154


Interpolatorische Quadraturformeln Definition: Newton-Cotes Formeln

Beispiele: Die Newton-Cotes Formeln mit


n = 1: Trapezregel

n = 2: Simpsonregel

n = 3: 38 -Regel
Z b      
b−a 2a + b a + 2b
f (x) dx ≈ f (a) + 3f + 3f + f (b)
a 8 3 3

Numerische Mathematik I 155


Interpolatorische Quadraturformeln Bemerkung:

4.1.4 Bemerkung: Ein einfacher Weg zur Berechnung der Gewichte: Anstatt der
Rb
Formel ωk = a Ln,k (x) dx kann man die Gewichte auch aus den
Exaktheitsforderungen bestimmen:
Rb
a 1 dx = ω0 + ω1 + ··· + ωn
Rb
a (x − a) dx = ω0 (x0 − a) + ω1 (x1 − a) + ··· + ωn (xn − a)
..
.
Rb
a
(x − a)n dx = ω0 (x0 − a)n + ω1 (x1 − a)n + ··· + ωn (xn − a)n

Dies ist ein lineares Gleichungssystem mit den Unbekannten ω0 , . . . , ωn und der
2 n+1
“rechten Seite” (b − a, (b−a)
2 , . . . , (b−a) T
n+1 ) , das für vorgegebene Knoten

x0 < x1 < · · · < xn

eindeutig lösbar ist (Transponierte der Vandermonde-Matrix!).

Numerische Mathematik I 156


Interpolatorische Quadraturformeln Bemerkung:

Bemerkung: Als Variante verwendet man auch die offenen Newton-Cotes Formeln,
bei denen die Randpunkte des Intervalls keine Knoten sind:
b−a
Knoten xk = a + k · , k = 0, 2, . . . , n,
n+2
Gewichte wie üblich (interpolatorisch)

Der Exaktheitsgrad der offenen Newton-Cotes Formeln mit n + 1 Knoten ist wieder

n für ungerades n, n + 1 für gerades n.

Als Beispiel für n = 0 ergibt sich die Mittelpunktsregel.

Numerische Mathematik I 157


Interpolatorische Quadraturformeln Satz: Fehler der Quadraturformel

Die Exaktheit für Polynome vom Grad ≤ n führt zur Darstellung des Fehlers der
Quadraturformel.

4.1.5 Satz: Fehler der Quadraturformel


Für eine interpolatorische Quadraturformel mit n + 1 Knoten x0 , . . . , xn gilt
Z b Z b n
Y
f (x) dx − In (f ; [a, b]) = f [x0 , . . . , xn , x] (x − xj ) dx.
a a j=0

Numerische Mathematik I 158


Interpolatorische Quadraturformeln Satz

4.1.6 Satz
Für eine (m + 1)-mal differenzierbare Funktion f : [a, b] → R gelten die folgenden
Aussagen:
a) n = 0, m = 1: Mittelpunktsregel
Z b  
a+b (b − a)3 ′′
f (x) dx − (b − a)f = f (ξ)
a 2 24

b) n = 1, m = 1: Trapezregel
Z b
b−a (b − a)3 ′′
f (x) dx − (f (a) + f (b)) = − f (ξ)
a 2 12

c) n = 2, m = 3: Simpsonregel
Z b    
b−a a+b (b − a)5 (4)
f (x) dx − f (a) + 4f + f (b) = − f (ξ)
a 6 2 2880

Hierbei ist m jeweils der Exaktheitsgrad der Quadraturformel, n + 1 ihre Knotenzahl und
ξ ∈ [a, b] ein geeigneter Punkt.

Numerische Mathematik I 159


Interpolatorische Quadraturformeln Satz

4.1.7 Satz
Für die summierten Quadraturformeln (Aufteilung von [a, b] in N Teilintervalle der Länge
h = (b − a)/N) gelten die folgenden Aussagen (mit geeignetem ξ ∈ [a, b]):
a’) summierte Mittelpunktsregel:
Z b
b−a (b − a)h2 ′′
f (x) dx − (f (y1 ) + f (y2 ) + · · · + f (yN )) = f (ξ)
a N 24

b’) summierte Trapezregel:


Z b
b−a
f (x) dx − (f (a) + 2f (a1 ) + 2f (a2 ) + · · · + 2f (aN−1 ) + f (b)) =
a 2N
(b − a)h2 ′′
− f (ξ)
12

c’) summierte Simpsonregel:


Z b b−a
f (x) dx − f (a) + 4f (y1 ) + 2f (a1 ) + 4f (y2 ) + · · ·
a 6N
!
(b − a)h4 (4)
+2f (aN−1 ) + 4f (yN ) + f (b) = − f (ξ)
2880

Numerische Mathematik I 160


Interpolatorische Quadraturformeln Beispiel:

4.1.8 Beispiel: Z
1 2
Zur Berechnung von dx = ln 2 = 0.693147... verwenden wir die Mittelpunkts-, Trapez- und
x 1
Simpsonregel sowie die summierten Regeln mit N = 2 und N = 4 Teilintervallen:

N 1 2 4

Trapezregel 0.75 0.708333 0.697024

Mittelpunktsregel 0.66 0.685714 0.691220

Simpsonregel 0.6944 0.693254 0.693154

Die Simpsonregel mit N = 4 Teilintervallen liefert schon 4 Nachkommastellen. Sie erfordert die
Auswertung des Integranden an 9 Stellen.

Numerische Mathematik I 161


Interpolatorische Quadraturformeln Beispiel:

In Anwendungen wird eine Fehlertoleranz ǫ > 0 vorgegeben. Dann muss die Anzahl N der
Teilintervalle so bestimmt werden, dass die summierte Quadraturformel die Fehlertoleranz nicht
überschreitet.
4.1.9 Beispiel: Wie groß muss N für die Aufteilung in Teilintervalle der Länge h = (b − a)/N
Z 2
1
gewählt sein, damit die Berechnung von dx = ln 2 = 0.693147... mit der Fehlertoleranz
1 x
ǫ = 10 −6 erfolgt:

1 1 2 6 24
, f ′ (x) = − 2 , f ′′ (x) = 3 , f ′′′ (x) = − 4 , f (4) (x) = 5 .
f (x) =
x x x x x
In den Fehlerdarstellungen verwenden wir
max |f ′′ (ξ)| = 2, max |f (4) (ξ)| = 24.
ξ∈[1,2] ξ∈[1,2]

(b − a)h2 ′′ 1 p
summierte Trapezregel: |f (ξ)| ≤ < 10−6 gilt für N > 106 /6, also
12 6N 2
N ≥ 409. Dafür sind N + 1 = 410 Auswertungen von f erforderlich.

(b − a)h2 ′′ 1 p
summierte Mittelpunktsregel: |f (ξ)| ≤ < 10−6 gilt für N > 106 /12,
24 12N 2
also N ≥ 289. Dafür sind N = 289 Auswertungen von f erforderlich.

(b − a)h4 (4) 1 p
summierte Simpsonregel: |f (ξ)| ≤ 4
< 10−6 gilt für N > 4 106 /120,
2880 120N
also N ≥ 10. Dafür sind 2N + 1 = 21 Auswertungen von f erforderlich.
Numerische Mathematik I 162
Gaußsche Quadraturformeln

4.2 Gaußsche Quadraturformeln

Frage: Können noch bessere Annäherungen erzielt werden, indem auch die Knoten
x0 , . . . , xn geschickter gewählt werden?
Als positive Antwort werden die Gauß-Formeln entwickelt. Dazu wählen wir die
Knoten xk und die Gewichte ωk , 0 ≤ k ≤ n, so, dass der größtmögliche
Exaktheitsgrad 2n + 1 erzielt wird.
Heuristische Überlegung: Die Exaktheitsbedingungen zum Grad 2n + 1 ergeben
2n + 2 nichtlineare Gleichungen mit den 2n + 2 Unbekannten x0 , . . . , xn und
ω0 , . . . , ωn .

Numerische Mathematik I 163


Gaußsche Quadraturformeln Beispiel:

4.2.1 Beispiel:: Zwei-punktige Gauß-Formel auf dem Intervall [−1, 1]


Die Exaktheitsbedingungen lauten
R1
−1
dx = 2 = ω0 + ω1
R1
−1
x dx = 0 = ω0 x 0 + ω1 x 1
R1 2
−1
x 2 dx = = ω0 x02 + ω1 x12
3
R1
−1 x 3 dx = 0 = ω0 x03 + ω1 x13

Die eindeutige Lösung lautet:


√ √
x0 = − 3/3, x1 = 3/3, ω0 = ω1 = 1,
also erhalten wir die Gauß-Quadraturformel
Z 1  √  √ 
f (x) dx ≈ f − 3/3 + f 3/3 .
−1

Beachte: Symmetrie der Knoten (x0 = −x1 ) und Gewichte (ω0 = ω1 )

Numerische Mathematik I 164


Gaußsche Quadraturformeln Satz: Nullstellen der Orthogonalpolynome

4.2.2 Satz: Nullstellen der Orthogonalpolynome


Das gewichtete Skalarprodukt h·, ·iω für C [−1, 1] erfülle die Voraussetzung der
3-Term Rekursion in Satz 3.3.12. Dann gilt für die Orthogonalpolynome e pn (mit
Höchstkoeffizient 1):
pn hat n einfache Nullstellen im Intervall (−1, 1),
e
(n)
−1 < ξ1 < · · · < ξn(n) < 1.

Die Nullstellen von e


pn−1 trennen die Nullstellen von e
pn ,
(n) (n−1) (n) (n−1) (n) (n−1)
−1 < ξ1 < ξ1 < ξ2 < ξ2 · · · < ξn−1 < ξn−1 < ξn(n) < 1.

Numerische Mathematik I 165


Gaußsche Quadraturformeln Satz: Nullstellen der Orthogonalpolynome

Beweis:
a) Setze N := {ξ ∈ (−1, 1) : ξ ist Nullstelle von e
pn mit ungerader Vielfachheit}.
Zu zeigen ist #N = n.
Q
Falls #N < n, setze q(x) = ξ∈N (x − ξ) ∈ Pn−1 . Dann hat qe pn keine Vorzeichenwechsel in
[−1, 1], also ist
Z 1
q(x)epn (x)ω(x) dx 6= 0,
−1
im Widerspruch zur Orthogonalität von e
pn zu Pn−1 .

Numerische Mathematik I 166


Gaußsche Quadraturformeln Satz: Nullstellen der Orthogonalpolynome

b) per Induktion nach n ∈ N0 :


Die Nullstellen von e
p0 (keine) trennen die von e
p1 .

Es gelte die Trennungseigenschaft für die Nullstellen von e pn−1 und epn . Insbesondere ist
(n)
pn−1 (ξj ) 6= 0 für j = 0, . . . , n. Die Drei-Term-Rekursion ergibt für e
e pn+1 ausgewertet an
den Nullstellen von e pn :
 
(n) (n) (n)
sign(epn+1 (ξj )) = sign − γn e pn−1 (ξj ) = −sign(e pn−1 (ξj )).
|{z}
>0

Wegen Höchstkoeffizient 1, Betrachtung der Grenzwerte für x → ±∞ und Lage aller


Nullstellen in (−1, 1) ist außerdem
pn+1 (1) = signe
signe pn−1 (1) = 1, pn−1 (−1) = (−1)n−1 .
pn+1 (−1) = signe
signe
Die Trennungseigenschaft der (einfachen) Nullstellen von e
pn−1 ergibt
(n)
pn−1 (ξj ) = (−1)n−j ,
signe j = 1, . . . , n,

also für e
pn+1
(n)
pn+1 (ξj ) = (−1)n+1−j ,
signe j = 1, . . . , n.
Daraus erhalten wir je eine Nullstelle von e
pn+1 in den offenen Intervallen
(n) (n) (n) (n) (n) (n)
(−1, ξ1 ), (ξ1 ), ξ2 ), . . . , (ξn−1 ), ξn ), (ξn ), 1).
Dies ist die Trennungseigenschaft der Nullstellen von e
pn und e
pn+1 .
Numerische Mathematik I 167
Gaußsche Quadraturformeln Satz: Gauß-Formeln auf [−1, 1]

4.2.3 Satz: Gauß-Formeln auf [−1, 1]


Für jedes n ∈ N0 existieren eindeutig bestimmte Knoten

−1 < x0 < x1 < · · · < xn < 1

und Gewichte ωk > 0, so dass die Quadraturformel


n
X
In (f ; [−1, 1]) = ωk f (xk )
k=0

den Exaktheitsgrad 2n + 1 hat.


Die Knoten sind die Nullstellen des Legendre-Polynoms Ln+1 in 3.3.11.
Die Gewichte erfüllen
Z 1 Z 1
ωk = Ln,k (x) dx = (Ln,k (x))2 dx > 0, k = 0, . . . , n.
−1 −1

Es gelten die Symmetrie-Eigenschaften

xk = −xn−k , ωk = ωn−k , k = 0, 1, . . . , n.
Numerische Mathematik I 168
Gaußsche Quadraturformeln Satz: Gauß-Formeln auf [−1, 1]

Beweis: Wir wählen x0 , . . . , xn als die Nullstellen von Ln+1 und die Gewichte für die
interpolatorische Q.f. als Z 1
ωk = Ln,k (x) dx.
−1

1. Diese Q.f. hat den Exaktheitsgrad 2n + 1:


Sei q ∈ P2n+1 . Polynomdivision ergibt
q = pLn+1 + r mit p, r ∈ Pn .
Damit ist einerseits
Z 1 Z 1 Z 1
q(x) dx = p(x)Ln+1 (x) dx + r (x) dx,
−1 −1 −1
| {z }
=0 wg.Orth.

und andererseits  
n
X
In (q) = ωk p(xk )Ln+1 (xk ) +r (xk ) = In (r ).
k=0
| {z }
=0 wg.Nullst.
Die Exaktheit von In für Pn folgt aus der Wahl der Gewichte, also insgesamt
Z 1 Z 1
q(x) dx = r (x) dx = In (r ) = In (q).
−1 −1

Numerische Mathematik I 169


Gaußsche Quadraturformeln Satz: Gauß-Formeln auf [−1, 1]

2. Darstellungen der Gewichte und Symmetrie:


Für die Gewichte dieser Q.f. folgt aus L2n,j ∈ P2n und der Exaktheit für P2n+1 sofort
n
X Z 1
ωk = ωj (Ln,j (xk ))2 = (Ln,j (x))2 dx,
j=0 −1

also insbesondere ωk > 0 für k = 0, . . . , n.


Symmetrie der Knoten: Ln+1 ist gerade bzw. ungerade.
Symmetrie der Gewichte: Ln,k (x) = Ln,n−k (−x)

3. Eindeutigkeit:
Sei In eine Q.f. mit paarweise verschiedenen Knoten x0 , . . . , xn , Gewichten
Q ω0 , . . . , ωn und
Exaktheitsgrad 2n + 1. Für das zugehörige Knotenpolynom w (x) = nj=0 (x − xj ) gilt
Z 1 n
X
w ∈ Pn+1 , q(x)w (x) dx = ωk q(xk ) w (xk ) = 0 für alle q ∈ Pn .
−1 | {z } k=0
| {z }
∈P2n+1 =0

Also ist w ⊥ Pn ; sein Höchstkoeffizient ist 1, und deshalb ist w = Ln+1 .

Numerische Mathematik I 170


Gaußsche Quadraturformeln Satz: Gauß-Formeln auf [−1, 1]

Beispiele:
L1 (x) = x hat die einzige Nullstelle x0 = 0. Die Ein-punktige Gauß-Formel ist
die Mittelpunktsregel.
√ √
L2 (x) = x 2 − 13 hat die Nullstellen x0 = − 3/3, x1 = 3/3. Dies sind die
Knoten der Zwei-punktigen Gauß-Formel in Bsp. 4.2.1.
p p
L3 (x) = x 3 − 35 x hat die Nullstellen x0 = − 3/5, x1 = 0 und x2 = 3/5.
Die Gewichte der 3-punktigen Gauß-Formel werden durch 3 lineare
Gleichungen aus den Exaktheitsbedingungen zum Grad 2 ermittelt:
5 8
ω0 = ω2 = , ω1 = .
9 9

Numerische Mathematik I 171


Gaußsche Quadraturformeln Satz: Fehlerformel

4.2.4 Satz: Fehlerformel


Die Gauß-Formel In , deren Knoten die Nullstellen des Legendre-Polynoms Ln+1
sind, erfüllt die Fehlerformel
Z 1
22n+3 ((n + 1)!)4 (2n+2)
f (x) dx − In (f ) = f (ξ), mit ξ ∈ [−1, 1]
−1 [(2n + 2)!]3 (2n + 3)

Numerische Mathematik I 172


Gaußsche Quadraturformeln Satz: Fehlerformel

Beweis: Wir betrachten neben der Gauß-Formel


n
X
In (f ; [−1, 1]) = ωk f (xk )
k=0

auch die Hermite-Quadraturformel


n
X
Jn (f ; [−1, 1]) = (αk f (xk ) + βk f ′ (xk )),
k=0

wobei die Wahl der Gewichte


Z 1 Z 1
αk = Hn,k (x) dx, βk = e n,k (x) dx
H
−1 −1

mit den Hermite-Grundpolynomen zu den doppelten Knoten x0 , . . . , xn erfolgt (siehe Übungsblatt


7). Die Quadraturformel Jn besteht aus Integration des Hermite-Interpolationspolynoms
n
X n
X
p(x) = f (xk )Hn,k (x) + f ′ (xk )H̃n,k (x)
k=0 k=0

und besitzt deshalb den Exaktheitsgrad 2n + 1. Die Integrale der Hermite-Grundpolynome sind
(wegen Exaktheitsgrad 2n + 1 sowohl von Jn als auch von In )
R
αk = Hn,k (x) dx = In (Hn,k ) = ωk ,
R−1,1
βk = e (x) dx = In (H
H e n,k ) = 0.
−1,1 n,k

Also stimmen In und Jn überein!


Numerische Mathematik I 173
Gaußsche Quadraturformeln Satz: Fehlerformel

Der Quadraturfehler wird mit Hilfe des Interpolationsfehlers der Hermite-Interpolation bestimmt
(vgl. Satz 3.1.15 mit dem erweiterten Knotenvektor)
Z b Z 1 n
Y
f (x) dx − In (f ) = f [x0 , x0 , . . . , xn , xn , x] (x − xj )2 dx.
a −1 j=0

Die dividierte Differenz hat die Ordnung 2n + 2, das Produkt ist (Ln+1 (x))2 . Der verallg. MWS
und die Normalisierung in Beispiel 3.3.11 ergeben
Z b Z 1
f (x) dx − In (f ) = f [x0 , x0 , . . . , xn , xn , ξ] (Ln+1 (x))2 dx
a −1
f (2n+2) (η) ((n + 1)!)4 22n+3
= .
(2n + 2)! ((2n + 2)!)2 2n + 3

Numerische Mathematik I 174


Gaußsche Quadraturformeln Satz: Fehlerformel

Bemerkung:
Z b
a) Die Gauß-Formel für f (x) dx erhält man durch Variablentransformation:
a

b−a
x ∈ [−1, 1] ⇐⇒ t(x) = a + (1 + x) ∈ [a, b]
2
ergibt die Knoten
b−a
tk = a + (1 + xk ) (xk Knoten in [−1, 1])
2
und die Gewichte
b−a
ηk = ωk .
2

b) Eine weitere Erhöhung der Genauigkeit erzielt man durch Aufteilen des
Intervalls in Teilintervalle (summierte Gauß-Formeln)

Numerische Mathematik I 175


Gaußsche Quadraturformeln Beispiel:

4.2.5 Beispiel:
Z 2
1
Für dx = ln 2 = 0.693147... verwenden wir die auf das Interall [1, 2]
1 x
transformierten Gauß-Formeln mit 2 bzw. 3 Knoten und die der summierten
Gauß-Formeln mit 2 Teilintervallen.
Zum Vergleich geben wir jeweils die Werte der Simpson-Regel an.

N 1 2

Simpsonregel 0.6944 0.693254

2-punkt. Gauß-Formel 0.692308 0.693077

3-punkt. Gauß-Formel 0.693122 0.693146

Numerische Mathematik I 176


Gaußsche Quadraturformeln Ergänzung: gewichtetes Integral

4.2.6 Ergänzung: gewichtetes Integral


Für das gewichtete Integral
Z 1
f (x) ω(x) dx
−1

erhält man den Exaktheitsgrad 2n + 1, wenn die Knoten xk als Nullstellen des
Orthogonalpolynoms e pn+1 vom Grad n + 1 und die Gewichte
Z 1 Z 1
ωk = Ln,k (x) ω(x) dx = (Ln,k (x))2 ω(x) dx > 0, k = 0, . . . , n,
−1 −1

gewählt werden. Dann ist der Quadraturfehler mit einem ξ ∈ (−1, 1)


Z 1 n
X Z 1
f (2n+2) (ξ)
f (x) ω(x) dx − ωk f (xk ) = pn+1 (x))2 ω(x) dx.
(e
−1 (2n + 2)! −1
k=0

Numerische Mathematik I 177


Gaußsche Quadraturformeln Ergänzung: gewichtetes Integral

Speziell für ω(x) = √1 :


1−x 2
Die Gauß-Tschebyscheff-Quadraturformeln haben die Knoten
 
π 2k + 1
xk = cos , k = 0, . . . , n,
2 n+1
dies sind die Nullstellen des Tschebyscheff-Polynoms Tn+1 erster Art (absteigend sortiert), siehe
Beispiel 3.3.9(b). Alle Gewichte
Z 1
dx
ωk = Ln,k (x) √
−1 1 − x2
π Pn
sind gleich, das liefert In (f ) = n+1 k=0 f (xk ).

Der Quadraturfehler hat die Darstellung


Z 1
dx π
f (x) √ − In (f ) = 2n+1 f (2n+2) (ξ), ξ ∈ (−1, 1).
−1 1 − x2 2 (2n + 2)!

Numerische Mathematik I 178


Gaußsche Quadraturformeln Ergänzung: gewichtetes Integral

Beweis: Wir zeigen die Exaktheit für Pn , indem wir den Quadraturfehler für Tm , m = 0, . . . , n,
berechnen. (
Z 1
dx π, m = 0,
Tm (x) √ =
−1 1 − x2 0, m ≥ 1
X n   (
π mπ 2k + 1 π, m = 0,
In (Tm ) = cos =
n + 1 k=0 2 n+1 0, 1 ≤ m ≤ n.
Die Summe 0 ergibt sich mit Hilfe der Eulerschen Formel aus
n   n  
mπ 2n+1 X mπ 2k + 1 X mπ 2n+2k+2 mπ 2n−2k
e i 2 n+1 cos = ei 2 n+1 + ei 2 n+1

k=0
2 n + 1 k=0

2n+1
X k 1 − e i 2mπ
= e imπ n+1 = mπ = 0.
k=0 1 − e i n+1
Die Konstante in der Fehlerdarstellung ergibt sich aus
Z 1 Z 1
dx 1 dx π
(w (x))2 √ = 2n (Tn+1 (x))2 √ = 2n+1 .
−1 1 − x2 2 −1 1 − x2 2
Man beachte hierbei, dass Tn+1 (x) = 2n w (x) gilt, weil Tn+1 den Höchstkoeffizienten 2n hat.

Numerische Mathematik I 179


Gaußsche Quadraturformeln Satz: Konvergenz der Gauß-Formeln

Die Positivität der Gewichte der Gauß-Formeln bewirkt, dass die Operator-Norm
von In auf C [−1, 1]

kIn k = sup{|In f | : f ∈ C [−1, 1], kf k∞ = 1}

unabhängig von n ∈ N0 den Wert


n
X Z 1
kIn k = ωk = ω(x) dx
k=0 −1

hat. Daraus resultiert eine wichtige Folgerung.

4.2.7 Satz: Konvergenz der Gauß-Formeln


Für die Gauß-Formel In mit n + 1 Knoten zum gewichteten Integral
R1
−1
f (x) ω(x) dx und jede Funktion f ∈ C [−1, 1] gilt
Z 1
lim In (f ) = f (x) ω(x) dx.
n→∞ −1

Numerische Mathematik I 180


Das Rombergsche Integrationsverfahren Hauptsatz: Entwicklung des Quadraturfehlers nach h-Potenzen

4.3 Das Rombergsche Integrationsverfahren

Die Verwendung der summierten Trapezregel zu verschiedenen Schrittweiten wird


als Ausgangspunkt zur Schrittweiten-Extrapolation nach Romberg verwendet.
4.3.1 Hauptsatz: Entwicklung des Quadraturfehlers nach h-Potenzen
Für f ∈ C 2m+2 [a, b] und N ∈ N, h = (b − a)/N, sei

X N−1
h h
a(h) = I1Σ (f ; N) = f (a) + h f (a + jh) + f (b).
2 2
j=1

Es gilt die Euler-Maclaurinsche Summenformel


Z b m
X B2k h2k
a(h) = f (x) dx + (f (2k−1) (b) − f (2k−1) (a)) +
a (2k)!
k=1

B2m+2 h2m+2 (2m+2)


(b − a) f (ξ)
(2m + 2)!

mit einem ξ ∈ [a, b] und den Bernoulli-Zahlen B2k .


Numerische Mathematik I 181
Das Rombergsche Integrationsverfahren Definition: Bernoulli-Polynome und Bernoulli-Zahlen

4.3.2 Definition: Bernoulli-Polynome und Bernoulli-Zahlen


Die Polynome bk ∈ Pk , k ∈ N0 , mit
Z 1
b0 (x) = 1, bk′ (x) = kbk−1 (x) und bk (x) dx = 0 für k ≥ 1,
0

heißen Bernoulli-Polynome, ihr Funktionswert Bk := bk (0) heißt Bernoulli-Zahl.

Bemerkung zu den Bernoulli-Zahlen:


a) Bernoulli-Zahlen sind
1 1 1 1
B0 = 1, B1 = − , B2 = , B3 = 0, B4 = − , B5 = 0, B6 =
2 6 30 42
b) Für k ∈ N ist B2k+1 = 0.

c) Rekursionsformel:
X
k−1
k Bj
B0 = 1, Bk = − für k = 1, 2, . . .
j=0
j k −j +1

d) Die Bernoulli-Zahlen treten vielfach auf, z.B. bei den Reihenwerten


X∞
1 (−1)k−1 π 2k 22k−1
2k
= B2k
n=1
n (2k)!

Numerische Mathematik I 182


Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:

4.3.3 Eigenschaften der Bernoulli-Polynome:


Z 1 Z 1
(i) Für k ≥ 2 ist bk (1) = bk (0) = Bk , weil bk′ (x) dx = k bk−1 (x) dx = 0 gilt.
0 0

(ii) Für k ∈ N0 gilt |B2k | = maxx∈[0,1] |b2k (x)|.


(ohne Beweis)

(iii) Für k ≥ 0 ist bk (x) = (−1)k bk (1 − x):

Beweis: klar für k = 0, per Induktion mit


d  
[bk − (−1)k bk (1 − ·)]′ (x) = k bk−1 (x) − (−1)k−1 bk−1 (1 − x) = 0,
dx
also bk − (−1)k bk (1 − ·) = ck mit einer Konstanten ck ∈ R für k ≥ 1. Man berechnet


 b (0) − bk (1) = 0 für gerades k ≥ 2 wegen (i),
 k
ck = Z 1 R


 (bk (x) + bk (1 − x)) dx = 0 für ungerades k ≥ 1 wegen 01 bk (x) dx = 0.
0

Numerische Mathematik I 183


Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:

(iv) Die Bernoulli-Polynome können mit Hilfe ihrer erzeugenden Funktion definiert werden:
X∞
te xt bk (x) k
= t , x ∈ R, 0 < |t| < 2π.
et − 1 k=0
k!

Mit dieser Darstellung lassen sich die meisten Eigenschaften der bk beweisen, indem man
einen Koeffizientenvergleich für Potenzreihen verwendet.
(Übungsaufgabe)

(v) Einsetzen von x = 0 ergibt die erzeugende Funktion der Bernoulli-Zahlen


X∞
t Bk k
= t , 0 < |t| < 2π.
et −1 k=0
k!

Numerische Mathematik I 184


Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:

Beweis von Satz 4.3.1: Für f ∈ C 2m+2 ([0, N]) folgt durch partielle Integration
Z 1 Z 1
f (x + ℓ) dx = b0 (x)f (x + ℓ) dx
0 0
Z 1
1
= (f (ℓ) + f (ℓ + 1)) − b1 (x)f ′ (x + ℓ) dx
2 0
Z 1
1 
= (f (ℓ) + f (ℓ + 1)) − B2 f ′ (ℓ + 1) − f ′ (ℓ) + b2 (x)f ′′ (x + ℓ) dx . . .
2 0

1 X B2k 
m+1 
= (f (ℓ) + f (ℓ + 1)) − f (2k−1) (ℓ + 1) − f (2k−1) (ℓ) +
2 k=1
(2k)!
Z 1
b2m+2 (x) (2m+2)
f (x + ℓ) dx.
0 (2m + 2)!

Hierbei wurde bk (1) = bk (0) = Bk und B2k+1 = 0 für k ≥ 1 benutzt.


Der letzte Summand zusammen mit dem letzten Integral ergibt
Z 1
b2m+2 (x) − B2m+2 (2m+2)
f (x + ℓ) dx.
0 (2m + 2)!
| {z }
festesVorzeichen

Das feste Vorzeichen von b2m+2 (x) − B2m+2 ergibt sich aus Eigenschaft (ii).

Numerische Mathematik I 185


Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:

Beim Zusammensetzen der Integrale von 0 bis N fallen einige Summanden wie bei einer
Teleskopsumme weg. Wir erhalten
Z
B2k  (2k−1) 
N N−1
X Xm
1
f (x) dx = (f (N) + f (0)) + f (j) − f (N) − f (2k−1) (0) +
0 2 j=1 k=1
(2k)!
Z N e
b2m+2 (x) − B2m+2 (2m+2)
f (x) dx.
0 (2m + 2)!

Dabei ist eb2m+2 die periodische Fortsetzung von b2m+2 vom Intervall [0, 1] auf [0, N], die
Differenz eb2m+2 (x) − B2m+2 hat also festes Vorzeichen auf [0, N]. Mit dem erweiterten
Mittelwertsatz folgt
Z N e Z N e
b2m+2 (x) − B2m+2 (2m+2) b2m+2 (x) − B2m+2
f (x) dx = f (2m+2) (ξ) dx
0 (2m + 2)! 0 (2m + 2)!
NB2m+2
= − f (2m+2) (ξ)
(2m + 2)!
mit ξ ∈ [0, N].
Die Aussage von Satz 4.3.1 folgt durch Koordinatentransformation von [0, N] auf [a, b].

Numerische Mathematik I 186


Das Rombergsche Integrationsverfahren Satz: Romberg-Verfahren

Direkte Anwendung des Extrapolationssatzes 3.2.3:


4.3.4 Satz: Romberg-Verfahren
Gegeben seien f ∈ C 2m+2 [a, b] sowie eine monoton wachsende Folge (Nj )j∈N0
natürlicher Zahlen mit
Nj
0< ≤ ρ < 1, j ∈ N0 .
Nj+1

Zu den Werten der summierten Trapezregel


b−a
aj,0 = a(hj ) = I1Σ (f ; Nj ), hj = ,
Nj

lauten die Einträge der Extrapolationstafel (nach Romberg)


aj,k−1 − aj−1,k−1
aj,k = aj,k−1 +  2 , k = 1, . . . , m, j = k, . . . , m.
hj −k
hj − 1

Numerische Mathematik I 187


Das Rombergsche Integrationsverfahren Satz: Romberg-Verfahren

Für 0 ≤ k ≤ m gilt
Z b
2k+2
f (x) dx − aj,k = O(hj−k ) für k ≤ j → ∞.
a

Speziell für hj = 2−j h0 gilt also


Z b
f (x) dx − aj,k = O(h02k+2 2−(2k+2)(j−k) ) für k ≤ j → ∞.
a

Numerische Mathematik I 188


Das Rombergsche Integrationsverfahren Folgerung für periodische Funktionen:

4.3.5 Folgerung für periodische Funktionen:


Für periodische Funktionen f ∈ C 2m+2 (R) mit Periodenlänge b − a fällt die summierte
Trapezregel zur Schrittweite h = (b − a)/N mit der Rechteckregel zusammen:

XN−1 X N−1
h h
I1Σ (f ; N) = f (a) + h f (a + jh) + f (b) = h f (a + jh).
2 j=1
2 j=0

Außerdem fallen die Terme f (2k−1) (b) − f (2k−1) (a) = 0 in der Euler-Maclaurinschen
Summenformel weg, d.h.
N−1
X Z b
a(h) = h f (a + jh) = f (x) dx + O(h2m+2 ).
j=0 a

Die Extrapolation führt hier zu keiner Verbesserung.


Ist sogar f ∈ C ∞ (R) periodisch mit Periodenlänge b − a, so konvergiert die Rechteckregel zur
R
Schrittweite h schneller gegen ab f (x) dx als jede Potenz von h. Die Verwendung komplizierterer
Quadraturformeln wäre hierfür eher schädlich!

Numerische Mathematik I 189


Praktische Aspekte der Integration

4.4 Praktische Aspekte der Integration

Ziel: Kriterien für die Wahl der Schrittweite der zusammengesetzten


Quadraturformeln zum Erzielen der Genauigkeit ǫ (oder eps)

Die a priori-Berechnung durch Verwendung der Fehlerformeln verlangt


Kenntnisse einer hohen Ableitung von f . Dies ist für praktische Zwecke nicht
sinnvoll.

Typ 1: Verwendung zweier Formeln zur Einschließung des Integralwerts.

Typ 2: Verwendung zweier Schrittweiten zur numerischen a


posteriori-Schätzung des Fehlers.

Numerische Mathematik I 190


Praktische Aspekte der Integration Typ 1: Einschließung des Integralwerts

4.4.1 Typ 1: Einschließung des Integralwerts


Die summierte Simpson-Regel zu N Unterteilungen, h = (b − a)/N, hat die Fehlerformel
Z b
(b − a)h4 (4)
f (x) dx − I2Σ (f ; N) = − f (ξ) mit ξ ∈ [a, b].
a 2880
Die summierte Form der offenen Newton-Cotes Formel (mit 3 Knoten pro Intervall) lautet
N−1
X
eI Σ (f ; N) = b − a (2f (ak + h/4) − f (ak + h/2) + 2f (ak + 3h/4)), ak = a + kh,
2
3N k=0

und hat die Fehlerformel


Z b
7(b − a)h4 (4)
f (x) dx − eI2Σ (f ; N) = f (ξ) mit ξ ∈ [a, b].
a 23040
Vergleich der beiden Fehlerterme zeigt, dass bei konstantem Vorzeichen von f (4) im Intervall
[a, b] der exakte Wert durch beide Quadraturformeln eingeschlossen wird, also
Z b Z b
eI Σ (f ; N) ≤ f (x) dx ≤ I2Σ (f ; N) oder I2Σ (f ; N) ≤ f (x) dx ≤ eI2Σ (f ; N)
2
a a

gilt.

Numerische Mathematik I 191


Praktische Aspekte der Integration Typ 1: Einschließung des Integralwerts

Die Voraussetzung an f (4)


ist i.A. nicht für [a, b] erfüllt, jedoch für (fast) alle Teilintervalle. Als
numerische Einschließung kann man daher
N−1
X Z b
min{I2 (f ; [ak , ak+1 ]), eI2 (f ; [ak , ak+1 ])} ≤ f (x) dx
k=0 a
N−1
X
≤ max{I2 (f ; [ak , ak+1 ]), eI2 (f ; [ak , ak+1 ])}
k=0

verwenden.

Numerische Mathematik I 192


Praktische Aspekte der Integration Typ 2: a-posteriori Fehlerschätzung

4.4.2 Typ 2: a-posteriori Fehlerschätzung


Es wird die Fehlerdarstellung der summierten Quadraturformel zur Schrittweite
h = (b − a)/N in der Form
Z b
f (x) dx − InΣ (f ; N) = cf hk + r (f , h)hk+1
a

mit einer Konstanten cf und einer beschränkten Funktion r (f , .) vorausgesetzt.

Es stehen zwei Werte InΣ (f ; N1 ) und InΣ (f ; N2 ) zu den Schrittweiten h1 > h2 > 0,
hj = (b − a)/Nj , j = 1, 2, zur Verfügung.
Dann kann die Konstante cf geschätzt werden durch Lösen von
 
cf (h1k − h2k ) = − InΣ (f ; N1 ) − InΣ (f ; N2 ) + O(h1k+1 ).

Für h1 ≪ 1 folgt
InΣ (f ; N1 ) − InΣ (f ; N2 )
cf ≈ − .
h1k − h2k
Die “richtige” Schrittweite h für die vorgegebene Toleranz ǫ ist nun
h = (ǫ/cf )1/k .

Numerische Mathematik I 193


Praktische Aspekte der Integration Praktischer Aspekt: adaptive Unterteilung

4.4.3 Praktischer Aspekt: adaptive Unterteilung


Bei den summierten Quadraturformeln zur Schrittweite h = (b − a)/N treten
in den Teilintervallen
p ganz unterschiedliche Fehlergrößen auf:
f (x) = |x − 0.7| ist bei x = 0.7 nicht differenzierbar, zu erwarten sind
große Fehler in diesem Abschnitt.

Daher wird die Schrittweite nicht überall halbiert, sondern nur in Intervallen
[ak , ak+1 ], wo ein Fehlerschätzer signalisiert, dass
Z ak+1
ǫ(ak+1 − ak )
f (x) dx − In (f , [ak , ak+1 ]) >
ak (b − a)

gilt. Dadurch entsteht eine adaptive Unterteilung von [a, b].

Numerische Mathematik I 194


Praktische Aspekte der Integration Beispiel:

4.4.4 Beispiel:
Der Matlab/Octave Befehl quad verwendet eine adaptive summierte Simpson-Regel. Die
Integration Z 1p
|x − 0.7| dx
0
zur Genauigkeit ǫ = 10−4 und Ausgabe der einzelnen Unterteilungsintervalle (mit Hilfe der
trace-Variablen) lautet

f=inline(’sqrt(abs(x-0.7))’);
trace=1;
q=quad(f,0,1,1e-4,trace)

Die alte Matlab-Version R12 verwendet adaptive Halbierung des Intervalls [0, 1] und erhält die
Zwischenpunkte
(a0 = 0), a1 = 0.5, a2 = 0.625, a3 = 0.6875, a4 = 0.75, (a5 = 1)
(→ Baumstruktur der Unterteilung). Die neue Version (Matlab R2012) verwendet eine
willkürliche Anfangs-Unterteilung in 3 Teilintervalle
[a, a + 0.27158(b − a)], [a + 0.27158(b − a), b − 0.27158(b − a)], [b − 0.27158(b − a), b]

und beginnt dann mit adaptiver Halbierung.

Numerische Mathematik I 195


Praktische Aspekte der Integration Beispiel:

adaptive summierte Simpson−Regel


0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Numerische Mathematik I 196

Das könnte Ihnen auch gefallen