Post on 19-Jul-2019
transcript
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: Naherungswerte fur bestimmte Integrale, wenn sich keine geschlossene Formder Stammfunktion finden lasst.
Beispiele:
Die Stammfunktion F (x) =
∫ x
0
e−t2 dt benotigt man in der
Wahrscheinlichkeitstheorie, um Werte der Normalverteilung zu berechnen.
Der Integralsinus
∫ x
0
sin t
tdt spielt eine Rolle in der Signalverarbeitung
Integral-Exponentialfunktion
∫ x
−∞
et
tdt und Integral-Logarithmus
∫ x
0
dt
ln t
Numerische Mathematik I 148
Interpolatorische Quadraturformeln
Moglichkeiten zur numerischen Berechnung von Naherungswerten:
Potenzreihe der Stammfunktion bilden und die Partialsummen an denIntegrationsgrenzen auswerten
Hier: Formeln fur Naherungsflachen verwenden (z.B. Trapezregel, KeplerscheFass-Regel)
Numerische Mathematik I 149
Interpolatorische Quadraturformeln Definition: Quadraturformel
4.1.1 Definition: Quadraturformel
Eine Naherung des bestimmten Integrals
∫ b
a
f (x) dx wird durch die
Quadraturformel
In(f ; [a, b]) =n∑
i=0
ωi f (xi )
angegeben, wobei
x0 < x1 < · · · < xn die Knoten (meistens in [a, b] gewahlt) und
ω0, ω1, . . . , ωn ∈ R die Gewichte
bezeichnen.Die Quadraturformel heißt interpolatorisch, wenn ihre Gewichte
ωk =
∫ b
a
Ln,k(x) dx , k = 0, . . . , n
sind, wobei Ln,k(x) =
n∏
j=0j 6=k
x − xj
xk − xjdie Lagrange-Grundpolynome bezeichnen.
Numerische Mathematik I 150
Interpolatorische Quadraturformeln Definition: Quadraturformel
Beachte: Fur jedes Polynom p ∈ Pn gilt
p(x) =n∑
k=0
p(xk )Ln,k (x). “p interpoliert sich selbst”
Also liefert jede interpolatorische Quadraturformel mit n + 1 Knoten
In(p; [a,b]) =n∑
k=0
ωkp(xk ) =
∫ b
a
p(x) dx ,
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), fallssie fur alle Polynome p ∈ Pm vom Grad kleiner oder gleich m den exakten Wert∫ b
ap(x) dx liefert.
Jede interpolatorische Quadraturformel (mit n + 1 Knoten) hat mindestens denExaktheitsgrad n und hochstens den Exaktheitsgrad 2n + 1.
Denn: mindestens Exaktheitsgrad n: Definition interpolatorisch
hochstens Exaktheitsgrad 2n + 1: fur q(x) =∏n
j=0(x − xj)2 gilt
∫ baq(x) dx > 0, aber In(q) = 0.
Numerische Mathematik I 151
Interpolatorische Quadraturformeln Beispiele:
4.1.2 Beispiele: Fur
∫ b
a
f (x) dx verwendete Quadraturformeln
a) (b − a)f
(a+ b
2
)Mittelpunktsregel
b)b − a
2(f (a) + f (b)) Trapezregel
c)b − a
6
(f (a) + 4f
(a+ b
2
)+ f (b)
)Simpsonregel
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
ak = a+ k · b − a
N, k = 0, 1, . . . ,N ,
und verwenden der Mittelpunkte yk = ak−1+ak2 , also
a = a0 < y1 < a1 < y2 < · · · < aN−1 < yN < aN = b
erhalt man die “summierten” Regeln:
a’) summierte Mittelpunktsregel
IΣ0 (f ;N) =b − a
N(f (y1) + f (y2) + · · ·+ f (yN ))
b’) summierte Trapezregel
IΣ1 (f ;N) = fracb − a2N (f (a) + 2f (a1) + 2f (a2) + · · ·+ 2f (aN−1) + f (b))
c’) summierte Simpsonregel
IΣ2 (f ;N) =b − a
6N(f (a) + 4f (y1) + 2f (a1) + 4f (y2) + · · ·+ 2f (aN−1) + 4f (yN) + f (b))
Numerische Mathematik I 153
Interpolatorische Quadraturformeln Definition: Newton-Cotes Formeln
Systematik zur Definition der “einfachen” Regeln: Gegeben sei das Integral∫ b
af (x) dx .
4.1.3 Definition: Newton-Cotes Formeln
Fur gegebenes n ∈ N wahlen wir
die Knoten xk = a+ k · b − a
n, k = 0, 1, . . . , n,
die Gewichte ωk =
∫ b
a
Ln,k(x) dx wie ublich fur eine interpolatorische
Quadraturformel.
Bemerkung: Der Exaktheitsgrad der Newton-Cotes Formeln mit n + 1 Knoten ist
n fur ungerades n, n + 1 fur 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
∫ b
a
f (x) dx ≈ b − a
8
(f (a) + 3f
(2a+ b
3
)+ 3f
(a+ 2b
3
)+ f (b)
)
Numerische Mathematik I 155
Interpolatorische Quadraturformeln Bemerkung:
4.1.4 Bemerkung: Ein einfacher Weg zur Berechnung der Gewichte: Anstatt der
Formel ωk =∫ b
aLn,k(x) dx kann man die Gewichte auch aus den
Exaktheitsforderungen bestimmen:
∫ b
a1 dx = ω0 + ω1 + · · · + ωn
∫ b
a(x − a) dx = ω0(x0 − a) + ω1(x1 − a) + · · · + ωn(xn − a)
...∫ b
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
“rechten Seite” (b − a, (b−a)2
2 , . . . , (b−a)n+1
n+1 )T , das fur vorgegebene Knoten
x0 < x1 < · · · < xn
eindeutig losbar 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:
Knoten xk = a+ k · b − a
n + 2, k = 0, 2, . . . , n,
Gewichte wie ublich (interpolatorisch)
Der Exaktheitsgrad der offenen Newton-Cotes Formeln mit n+1 Knoten ist wieder
n fur ungerades n, n + 1 fur gerades n.
Als Beispiel fur n = 0 ergibt sich die Mittelpunktsregel.
Numerische Mathematik I 157
Interpolatorische Quadraturformeln Satz: Fehler der Quadraturformel
Die Exaktheit fur Polynome vom Grad ≤ n fuhrt zur Darstellung des Fehlers derQuadraturformel.
4.1.5 Satz: Fehler der Quadraturformel
Fur eine interpolatorische Quadraturformel mit n + 1 Knoten x0, . . . , xn gilt
∫ b
a
f (x) dx − In(f ; [a, b]) =
∫ b
a
f [x0, . . . , xn, x ]n∏
j=0
(x − xj) dx .
Numerische Mathematik I 158
Interpolatorische Quadraturformeln Satz
4.1.6 Satz
Fur eine (m + 1)-mal differenzierbare Funktion f : [a, b] → R gelten die folgendenAussagen:
a) n = 0, m = 1: Mittelpunktsregel
∫ b
a
f (x) dx − (b − a)f
(a+ b
2
)=
(b − a)3
24f ′′(ξ)
b) n = 1, m = 1: Trapezregel
∫ b
a
f (x) dx − b − a
2(f (a) + f (b)) = − (b − a)3
12f ′′(ξ)
c) n = 2, m = 3: Simpsonregel
∫ b
a
f (x) dx − b − a
6
(f (a) + 4f
(a+ b
2
)+ f (b)
)= − (b − a)5
2880f (4)(ξ)
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
Fur die summierten Quadraturformeln (Aufteilung von [a,b] in N Teilintervalle der Langeh = (b − a)/N) gelten die folgenden Aussagen (mit geeignetem ξ ∈ [a,b]):
a’) summierte Mittelpunktsregel:
∫ b
a
f (x) dx − b − a
N(f (y1) + f (y2) + · · ·+ f (yN )) =
(b − a)h2
24f ′′(ξ)
b’) summierte Trapezregel:
∫ b
a
f (x) dx − b − a
2N(f (a) + 2f (a1) + 2f (a2) + · · ·+ 2f (aN−1) + f (b)) =
− (b − a)h2
12f ′′(ξ)
c’) summierte Simpsonregel:
∫ b
a
f (x) dx − b − a
6N
(f (a) + 4f (y1) + 2f (a1) + 4f (y2) + · · ·
+2f (aN−1) + 4f (yN) + f (b)
)= − (b − a)h4
2880f (4)(ξ)
Numerische Mathematik I 160
Interpolatorische Quadraturformeln Beispiel:
4.1.8 Beispiel:
Zur Berechnung von
∫ 2
1
1
xdx = ln 2 = 0.693147... verwenden wir die Mittelpunkts-, Trapez- und
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 dieAuswertung des Integranden an 9 Stellen.
Numerische Mathematik I 161
Interpolatorische Quadraturformeln Beispiel:
In Anwendungen wird eine Fehlertoleranz ǫ > 0 vorgegeben. Dann muss die Anzahl N derTeilintervalle so bestimmt werden, dass die summierte Quadraturformel die Fehlertoleranz nichtuberschreitet.
4.1.9 Beispiel: Wie groß muss N fur die Aufteilung in Teilintervalle der Lange h = (b − a)/N
gewahlt sein, damit die Berechnung von
∫ 2
1
1
xdx = ln 2 = 0.693147... mit der Fehlertoleranz
ǫ = 10−6 erfolgt:
f (x) =1
x, f ′(x) = − 1
x2, f ′′(x) =
2
x3, f ′′′(x) = − 6
x4, f (4)(x) =
24
x5.
In den Fehlerdarstellungen verwenden wir
maxξ∈[1,2]
|f ′′(ξ)| = 2, maxξ∈[1,2]
|f (4)(ξ)| = 24.
summierte Trapezregel:(b − a)h2
12|f ′′(ξ)| ≤ 1
6N2< 10−6 gilt fur N >
√106/6, also
N ≥ 409. Dafur sind N + 1 = 410 Auswertungen von f erforderlich.
summierte Mittelpunktsregel:(b − a)h2
24|f ′′(ξ)| ≤ 1
12N2< 10−6 gilt fur N >
√106/12,
also N ≥ 289. Dafur sind N = 289 Auswertungen von f erforderlich.
summierte Simpsonregel:(b − a)h4
2880|f (4)(ξ)| ≤ 1
120N4< 10−6 gilt fur N > 4
√106/120,
also N ≥ 10. Dafur sind 2N + 1 = 21 Auswertungen von f erforderlich.
Numerische Mathematik I 162
Gaußsche Quadraturformeln
4.2 Gaußsche Quadraturformeln
Frage: Konnen noch bessere Annaherungen erzielt werden, indem auch die Knotenx0, . . . , xn geschickter gewahlt werden?
Als positive Antwort werden die Gauß-Formeln entwickelt. Dazu wahlen wir dieKnoten xk und die Gewichte ωk , 0 ≤ k ≤ n, so, dass der großtmoglicheExaktheitsgrad 2n+ 1 erzielt wird.
Heuristische Uberlegung: Die Exaktheitsbedingungen zum Grad 2n + 1 ergeben2n+ 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∫ 1−1
dx = 2 = ω0 + ω1
∫ 1−1
x dx = 0 = ω0x0 + ω1x1
∫ 1−1
x2 dx =2
3= ω0x
20 + ω1x
21
∫ 1−1 x
3 dx = 0 = ω0x30 + ω1x
31
Die eindeutige Losung lautet:
x0 = −√3/3, x1 =
√3/3, ω0 = ω1 = 1,
also erhalten wir die Gauß-Quadraturformel∫ 1
−1f (x) dx ≈ f
(−√3/3)+ f
(√3/3).
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 〈·, ·〉ω fur C [−1, 1] erfulle die Voraussetzung der3-Term Rekursion in Satz 3.3.12. Dann gilt fur die Orthogonalpolynome pn (mitHochstkoeffizient 1):
pn hat n einfache Nullstellen im Intervall (−1, 1),
−1 < ξ(n)1 < · · · < ξ(n)n < 1.
Die Nullstellen von pn−1 trennen die Nullstellen von pn,
−1 < ξ(n)1 < ξ
(n−1)1 < ξ
(n)2 < ξ
(n−1)2 · · · < ξ
(n)n−1 < ξ
(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 pn mit ungerader Vielfachheit}.Zu zeigen ist #N = n.
Falls #N < n, setze q(x) =∏
ξ∈N (x − ξ) ∈ Pn−1. Dann hat qpn keine Vorzeichenwechsel in
[−1, 1], also ist ∫ 1
−1q(x)pn(x)ω(x) dx 6= 0,
im Widerspruch zur Orthogonalitat von pn zu Pn−1.
Numerische Mathematik I 166
Gaußsche Quadraturformeln Satz: Nullstellen der Orthogonalpolynome
b) per Induktion nach n ∈ N0 :
Die Nullstellen von p0 (keine) trennen die von p1.
Es gelte die Trennungseigenschaft fur die Nullstellen von pn−1 und pn. Insbesondere ist
pn−1(ξ(n)j
) 6= 0 fur j = 0, . . . , n. Die Drei-Term-Rekursion ergibt fur pn+1 ausgewertet an
den Nullstellen von pn:
sign(pn+1(ξ(n)j
)) = sign
(− γn︸︷︷︸
>0
pn−1(ξ(n)j
)
)= −sign(pn−1(ξ
(n)j
)).
Wegen Hochstkoeffizient 1, Betrachtung der Grenzwerte fur x → ±∞ und Lage allerNullstellen in (−1, 1) ist außerdem
signpn+1(1) = signpn−1(1) = 1, signpn+1(−1) = signpn−1(−1) = (−1)n−1.
Die Trennungseigenschaft der (einfachen) Nullstellen von pn−1 ergibt
signpn−1(ξ(n)j
) = (−1)n−j , j = 1, . . . , n,
also fur pn+1
signpn+1(ξ(n)j
) = (−1)n+1−j , j = 1, . . . , n.
Daraus erhalten wir je eine Nullstelle von pn+1 in den offenen Intervallen
(−1, ξ(n)1 ), (ξ
(n)1 ), ξ
(n)2 ), . . . , (ξ
(n)n−1), ξ
(n)n ), (ξ
(n)n ), 1).
Dies ist die Trennungseigenschaft der Nullstellen von pn und pn+1.
Numerische Mathematik I 167
Gaußsche Quadraturformeln Satz: Gauß-Formeln auf [−1, 1]
4.2.3 Satz: Gauß-Formeln auf [−1, 1]
Fur jedes n ∈ N0 existieren eindeutig bestimmte Knoten
−1 < x0 < x1 < · · · < xn < 1
und Gewichte ωk > 0, so dass die Quadraturformel
In(f ; [−1, 1]) =n∑
k=0
ωk f (xk )
den Exaktheitsgrad 2n + 1 hat.
Die Knoten sind die Nullstellen des Legendre-Polynoms Ln+1 in 3.3.11.
Die Gewichte erfullen
ωk =
∫ 1
−1
Ln,k(x) dx =
∫ 1
−1
(Ln,k(x))2 dx > 0, k = 0, . . . , n.
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 wahlen x0, . . . , xn als die Nullstellen von Ln+1 und die Gewichte fur dieinterpolatorische Q.f. als
ωk =
∫ 1
−1Ln,k(x) dx .
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∫ 1
−1q(x) dx =
∫ 1
−1p(x)Ln+1(x) dx
︸ ︷︷ ︸=0 wg.Orth.
+
∫ 1
−1r(x) dx ,
und andererseits
In(q) =
n∑
k=0
ωk
(p(xk)Ln+1(xk)︸ ︷︷ ︸=0 wg.Nullst.
+r(xk)
)= In(r).
Die Exaktheit von In fur Pn folgt aus der Wahl der Gewichte, also insgesamt∫ 1
−1q(x) dx =
∫ 1
−1r(x) dx = In(r) = In(q).
Numerische Mathematik I 169
Gaußsche Quadraturformeln Satz: Gauß-Formeln auf [−1, 1]
2. Darstellungen der Gewichte und Symmetrie:
Fur die Gewichte dieser Q.f. folgt aus L2n,j ∈ P2n und der Exaktheit fur P2n+1 sofort
ωk =n∑
j=0
ωj (Ln,j (xk))2 =
∫ 1
−1(Ln,j (x))
2 dx ,
also insbesondere ωk > 0 fur 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 ω0, . . . , ωn undExaktheitsgrad 2n + 1. Fur das zugehorige Knotenpolynom w(x) =
∏nj=0(x − xj ) gilt
w ∈ Pn+1,
∫ 1
−1q(x)w(x)︸ ︷︷ ︸∈P2n+1
dx =n∑
k=0
ωkq(xk)w(xk)︸ ︷︷ ︸=0
= 0 fur alle q ∈ Pn.
Also ist w ⊥ Pn; sein Hochstkoeffizient 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 istdie Mittelpunktsregel.
L2(x) = x2 − 13 hat die Nullstellen x0 = −
√3/3, x1 =
√3/3. Dies sind die
Knoten der Zwei-punktigen Gauß-Formel in Bsp. 4.2.1.
L3(x) = x3 − 35 x hat die Nullstellen x0 = −
√3/5, x1 = 0 und x2 =
√3/5.
Die Gewichte der 3-punktigen Gauß-Formel werden durch 3 lineareGleichungen aus den Exaktheitsbedingungen zum Grad 2 ermittelt:
ω0 = ω2 =5
9, ω1 =
8
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, erfullt die Fehlerformel
∫ 1
−1
f (x) dx − In(f ) =22n+3((n + 1)!)4
[(2n+ 2)!]3(2n + 3)f (2n+2)(ξ), mit ξ ∈ [−1, 1]
Numerische Mathematik I 172
Gaußsche Quadraturformeln Satz: Fehlerformel
Beweis: Wir betrachten neben der Gauß-Formel
In(f ; [−1, 1]) =n∑
k=0
ωk f (xk)
auch die Hermite-Quadraturformel
Jn(f ; [−1, 1]) =n∑
k=0
(αk f (xk) + βk f′(xk)),
wobei die Wahl der Gewichte
αk =
∫ 1
−1Hn,k(x) dx , βk =
∫ 1
−1Hn,k(x) dx
mit den Hermite-Grundpolynomen zu den doppelten Knoten x0, . . . , xn erfolgt (siehe Ubungsblatt7). Die Quadraturformel Jn besteht aus Integration des Hermite-Interpolationspolynoms
p(x) =n∑
k=0
f (xk)Hn,k(x) +n∑
k=0
f ′(xk )Hn,k(x)
und besitzt deshalb den Exaktheitsgrad 2n + 1. Die Integrale der Hermite-Grundpolynome sind(wegen Exaktheitsgrad 2n + 1 sowohl von Jn als auch von In)
αk =∫−1,1 Hn,k(x) dx = In(Hn,k ) = ωk ,
βk =∫−1,1
Hn,k(x) dx = In(Hn,k ) = 0.
Also stimmen In und Jn uberein!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)
∫ b
a
f (x) dx − In(f ) =
∫ 1
−1f [x0, x0, . . . , xn, xn, x ]
n∏
j=0
(x − xj)2 dx .
Die dividierte Differenz hat die Ordnung 2n + 2, das Produkt ist (Ln+1(x))2. Der verallg. MWSund die Normalisierung in Beispiel 3.3.11 ergeben
∫ b
a
f (x) dx − In(f ) = f [x0, x0, . . . , xn, xn, ξ]
∫ 1
−1(Ln+1(x))
2 dx
=f (2n+2)(η)
(2n + 2)!
((n + 1)!)4
((2n + 2)!)222n+3
2n + 3.
Numerische Mathematik I 174
Gaußsche Quadraturformeln Satz: Fehlerformel
Bemerkung:
a) Die Gauß-Formel fur
∫ b
a
f (x) dx erhalt man durch Variablentransformation:
x ∈ [−1, 1] ⇐⇒ t(x) = a+b − a
2(1 + x) ∈ [a, b]
ergibt die Knoten
tk = a +b − a
2(1 + xk ) (xk Knoten in [−1, 1])
und die Gewichte
ηk =b − a
2ωk .
b) Eine weitere Erhohung der Genauigkeit erzielt man durch Aufteilen desIntervalls in Teilintervalle (summierte Gauß-Formeln)
Numerische Mathematik I 175
Gaußsche Quadraturformeln Beispiel:
4.2.5 Beispiel:
Fur
∫ 2
1
1
xdx = ln 2 = 0.693147... verwenden wir die auf das Interall [1, 2]
transformierten Gauß-Formeln mit 2 bzw. 3 Knoten und die der summiertenGauß-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 Erganzung: gewichtetes Integral
4.2.6 Erganzung: gewichtetes IntegralFur das gewichtete Integral ∫ 1
−1
f (x)ω(x) dx
erhalt man den Exaktheitsgrad 2n + 1, wenn die Knoten xk als Nullstellen desOrthogonalpolynoms pn+1 vom Grad n + 1 und die Gewichte
ωk =
∫ 1
−1
Ln,k(x)ω(x) dx =
∫ 1
−1
(Ln,k(x))2 ω(x) dx > 0, k = 0, . . . , n,
gewahlt werden. Dann ist der Quadraturfehler mit einem ξ ∈ (−1, 1)
∫ 1
−1
f (x)ω(x) dx −n∑
k=0
ωk f (xk ) =f (2n+2)(ξ)
(2n + 2)!
∫ 1
−1
(pn+1(x))2 ω(x) dx .
Numerische Mathematik I 177
Gaußsche Quadraturformeln Erganzung: gewichtetes Integral
Speziell fur ω(x) = 1√1−x2
:
Die Gauß-Tschebyscheff-Quadraturformeln haben die Knoten
xk = cos
(π
2
2k + 1
n + 1
), k = 0, . . . , n,
dies sind die Nullstellen des Tschebyscheff-Polynoms Tn+1 erster Art (absteigend sortiert), sieheBeispiel 3.3.9(b). Alle Gewichte
ωk =
∫ 1
−1Ln,k (x)
dx√1− x2
sind gleich, das liefert In(f ) =π
n+1
∑nk=0 f (xk ).
Der Quadraturfehler hat die Darstellung∫ 1
−1f (x)
dx√1− x2
− In(f ) =π
22n+1(2n + 2)!f (2n+2)(ξ), ξ ∈ (−1, 1).
Numerische Mathematik I 178
Gaußsche Quadraturformeln Erganzung: gewichtetes Integral
Beweis: Wir zeigen die Exaktheit fur Pn, indem wir den Quadraturfehler fur Tm , m = 0, . . . , n,berechnen. ∫ 1
−1Tm(x)
dx√1− x2
=
{π, m = 0,
0, m ≥ 1
In(Tm) =π
n + 1
n∑
k=0
cos
(mπ
2
2k + 1
n + 1
)=
{π, m = 0,
0, 1 ≤ m ≤ n.
Die Summe 0 ergibt sich mit Hilfe der Eulerschen Formel aus
ei mπ
22n+1n+1
n∑
k=0
cos
(mπ
2
2k + 1
n + 1
)=
n∑
k=0
(ei mπ
22n+2k+2
n+1 + ei mπ
22n−2kn+1
)
=2n+1∑
k=0
eimπ k
n+1 =1− e i2mπ
1− ei mπ
n+1
= 0.
Die Konstante in der Fehlerdarstellung ergibt sich aus∫ 1
−1(w(x))2
dx√1− x2
=1
22n
∫ 1
−1(Tn+1(x))
2 dx√1− x2
=π
22n+1.
Man beachte hierbei, dass Tn+1(x) = 2nw(x) gilt, weil Tn+1 den Hochstkoeffizienten 2n hat.
Numerische Mathematik I 179
Gaußsche Quadraturformeln Satz: Konvergenz der Gauß-Formeln
Die Positivitat der Gewichte der Gauß-Formeln bewirkt, dass die Operator-Normvon In auf C [−1, 1]
‖In‖ = sup{|Inf | : f ∈ C [−1, 1], ‖f ‖∞ = 1}
unabhangig von n ∈ N0 den Wert
‖In‖ =
n∑
k=0
ωk =
∫ 1
−1
ω(x) dx
hat. Daraus resultiert eine wichtige Folgerung.
4.2.7 Satz: Konvergenz der Gauß-Formeln
Fur die Gauß-Formel In mit n + 1 Knoten zum gewichteten Integral∫ 1
−1f (x)ω(x) dx und jede Funktion f ∈ C [−1, 1] gilt
limn→∞
In(f ) =
∫ 1
−1
f (x)ω(x) dx .
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 wirdals Ausgangspunkt zur Schrittweiten-Extrapolation nach Romberg verwendet.
4.3.1 Hauptsatz: Entwicklung des Quadraturfehlers nach h-Potenzen
Fur f ∈ C 2m+2[a, b] und N ∈ N, h = (b − a)/N , sei
a(h) = IΣ1 (f ;N) =h
2f (a) + h
N−1∑
j=1
f (a+ jh) +h
2f (b).
Es gilt die Euler-Maclaurinsche Summenformel
a(h) =
∫ b
a
f (x) dx +
m∑
k=1
B2kh2k
(2k)!(f (2k−1)(b)− f (2k−1)(a)) +
(b − a)B2m+2h
2m+2
(2m+ 2)!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
b0(x) = 1, b′k(x) = kbk−1(x) und
∫ 1
0bk (x) dx = 0 fur k ≥ 1,
heißen Bernoulli-Polynome, ihr Funktionswert Bk := bk(0) heißt Bernoulli-Zahl.
Bemerkung zu den Bernoulli-Zahlen:
a) Bernoulli-Zahlen sind
B0 = 1, B1 = −1
2, B2 =
1
6, B3 = 0, B4 = − 1
30, B5 = 0, B6 =
1
42
b) Fur k ∈ N ist B2k+1 = 0.
c) Rekursionsformel:
B0 = 1, Bk = −k−1∑
j=0
(kj
) Bj
k − j + 1fur k = 1, 2, . . .
d) Die Bernoulli-Zahlen treten vielfach auf, z.B. bei den Reihenwerten
∞∑
n=1
1
n2k=
(−1)k−1π2k22k−1
(2k)!B2k
Numerische Mathematik I 182
Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:
4.3.3 Eigenschaften der Bernoulli-Polynome:
(i) Fur k ≥ 2 ist bk (1) = bk (0) = Bk , weil
∫ 1
0b′k (x) dx = k
∫ 1
0bk−1(x) dx = 0 gilt.
(ii) Fur k ∈ N0 gilt |B2k | = maxx∈[0,1] |b2k(x)|.(ohne Beweis)
(iii) Fur k ≥ 0 ist bk (x) = (−1)kbk(1 − x):
Beweis: klar fur k = 0, per Induktion mit
d
dx[bk − (−1)kbk (1− ·)]′(x) = k
(bk−1(x)− (−1)k−1bk−1(1− x)
)= 0,
also bk − (−1)kbk(1 − ·) = ck mit einer Konstanten ck ∈ R fur k ≥ 1. Man berechnet
ck =
bk(0) − bk(1) = 0 fur gerades k ≥ 2 wegen (i),
∫ 1
0(bk (x) + bk(1 − x)) dx = 0 fur ungerades k ≥ 1 wegen
∫ 10bk (x) dx = 0.
Numerische Mathematik I 183
Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:
(iv) Die Bernoulli-Polynome konnen mit Hilfe ihrer erzeugenden Funktion definiert werden:
text
et − 1=
∞∑
k=0
bk (x)
k!tk , x ∈ R, 0 < |t| < 2π.
Mit dieser Darstellung lassen sich die meisten Eigenschaften der bk beweisen, indem maneinen Koeffizientenvergleich fur Potenzreihen verwendet.
(Ubungsaufgabe)
(v) Einsetzen von x = 0 ergibt die erzeugende Funktion der Bernoulli-Zahlen
t
et − 1=
∞∑
k=0
Bk
k!tk , 0 < |t| < 2π.
Numerische Mathematik I 184
Das Rombergsche Integrationsverfahren Eigenschaften der Bernoulli-Polynome:
Beweis von Satz 4.3.1: Fur f ∈ C2m+2([0,N]) folgt durch partielle Integration
∫ 1
0f (x + ℓ) dx =
∫ 1
0b0(x)f (x + ℓ) dx
=1
2(f (ℓ) + f (ℓ+ 1)) −
∫ 1
0b1(x)f
′(x + ℓ) dx
=1
2(f (ℓ) + f (ℓ+ 1)) − B2
(f ′(ℓ+ 1)− f ′(ℓ)
)+
∫ 1
0b2(x)f
′′(x + ℓ) dx . . .
=1
2(f (ℓ) + f (ℓ+ 1)) −
m+1∑
k=1
B2k
(2k)!
(f (2k−1)(ℓ+ 1) − f (2k−1)(ℓ)
)+
∫ 1
0
b2m+2(x)
(2m + 2)!f (2m+2)(x + ℓ) dx .
Hierbei wurde bk (1) = bk (0) = Bk und B2k+1 = 0 fur k ≥ 1 benutzt.
Der letzte Summand zusammen mit dem letzten Integral ergibt∫ 1
0
b2m+2(x)− B2m+2
(2m + 2)!︸ ︷︷ ︸
festesVorzeichen
f (2m+2)(x + ℓ) dx .
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 einerTeleskopsumme weg. Wir erhalten
∫ N
0f (x) dx =
1
2(f (N) + f (0)) +
N−1∑
j=1
f (j)−m∑
k=1
B2k
(2k)!
(f (2k−1)(N) − f (2k−1)(0)
)+
∫ N
0
b2m+2(x)− B2m+2
(2m + 2)!f (2m+2)(x) dx .
Dabei ist b2m+2 die periodische Fortsetzung von b2m+2 vom Intervall [0, 1] auf [0,N], die
Differenz b2m+2(x)− B2m+2 hat also festes Vorzeichen auf [0,N]. Mit dem erweitertenMittelwertsatz folgt
∫ N
0
b2m+2(x) − B2m+2
(2m + 2)!f (2m+2)(x) dx = f (2m+2)(ξ)
∫ N
0
b2m+2(x)− B2m+2
(2m + 2)!dx
= − NB2m+2
(2m + 2)!f (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
naturlicher Zahlen mit
0 <Nj
Nj+1≤ ρ < 1, j ∈ N0.
Zu den Werten der summierten Trapezregel
aj,0 = a(hj) = IΣ1 (f ;Nj), hj =b − a
Nj
,
lauten die Eintrage der Extrapolationstafel (nach Romberg)
aj,k = aj,k−1 +aj,k−1 − aj−1,k−1(
hj−k
hj
)2
− 1
, k = 1, . . . ,m, j = k , . . . ,m.
Numerische Mathematik I 187
Das Rombergsche Integrationsverfahren Satz: Romberg-Verfahren
Fur 0 ≤ k ≤ m gilt
∫ b
a
f (x) dx − aj,k = O(h2k+2j−k ) fur k ≤ j → ∞.
Speziell fur hj = 2−jh0 gilt also
∫ b
a
f (x) dx − aj,k = O(h2k+20 2−(2k+2)(j−k)) fur k ≤ j → ∞.
Numerische Mathematik I 188
Das Rombergsche Integrationsverfahren Folgerung fur periodische Funktionen:
4.3.5 Folgerung fur periodische Funktionen:Fur periodische Funktionen f ∈ C2m+2(R) mit Periodenlange b − a fallt die summierteTrapezregel zur Schrittweite h = (b − a)/N mit der Rechteckregel zusammen:
IΣ1 (f ;N) =h
2f (a) + h
N−1∑
j=1
f (a + jh) +h
2f (b) = h
N−1∑
j=0
f (a+ jh).
Außerdem fallen die Terme f (2k−1)(b) − f (2k−1)(a) = 0 in der Euler-MaclaurinschenSummenformel weg, d.h.
a(h) = h
N−1∑
j=0
f (a + jh) =
∫ b
a
f (x) dx +O(h2m+2).
Die Extrapolation fuhrt hier zu keiner Verbesserung.
Ist sogar f ∈ C∞(R) periodisch mit Periodenlange b − a, so konvergiert die Rechteckregel zur
Schrittweite h schneller gegen∫ b
af (x) dx als jede Potenz von h. Die Verwendung komplizierterer
Quadraturformeln ware hierfur eher schadlich!
Numerische Mathematik I 189
Praktische Aspekte der Integration
4.4 Praktische Aspekte der Integration
Ziel: Kriterien fur die Wahl der Schrittweite der zusammengesetztenQuadraturformeln zum Erzielen der Genauigkeit ǫ (oder eps)
Die a priori-Berechnung durch Verwendung der Fehlerformeln verlangtKenntnisse einer hohen Ableitung von f . Dies ist fur praktische Zwecke nichtsinnvoll.
Typ 1: Verwendung zweier Formeln zur Einschließung des Integralwerts.
Typ 2: Verwendung zweier Schrittweiten zur numerischen a
posteriori-Schatzung des Fehlers.
Numerische Mathematik I 190
Praktische Aspekte der Integration Typ 1: Einschließung des Integralwerts
4.4.1 Typ 1: Einschließung des IntegralwertsDie summierte Simpson-Regel zu N Unterteilungen, h = (b − a)/N, hat die Fehlerformel
∫ b
a
f (x) dx − IΣ2 (f ;N) = − (b − a)h4
2880f (4)(ξ) mit ξ ∈ [a,b].
Die summierte Form der offenen Newton-Cotes Formel (mit 3 Knoten pro Intervall) lautet
IΣ2 (f ;N) =b − a
3N
N−1∑
k=0
(2f (ak + h/4) − f (ak + h/2) + 2f (ak + 3h/4)), ak = a+ kh,
und hat die Fehlerformel∫ b
a
f (x) dx − IΣ2 (f ;N) =7(b − a)h4
23040f (4)(ξ) mit ξ ∈ [a, b].
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
IΣ2 (f ;N) ≤∫ b
a
f (x) dx ≤ IΣ2 (f ;N) oder IΣ2 (f ;N) ≤∫ b
a
f (x) dx ≤ IΣ2 (f ;N)
gilt.
Numerische Mathematik I 191
Praktische Aspekte der Integration Typ 1: Einschließung des Integralwerts
Die Voraussetzung an f (4) ist i.A. nicht fur [a,b] erfullt, jedoch fur (fast) alle Teilintervalle. Alsnumerische Einschließung kann man daher
N−1∑
k=0
min{I2(f ; [ak , ak+1]), I2(f ; [ak , ak+1])} ≤∫ b
a
f (x) dx
≤N−1∑
k=0
max{I2(f ; [ak , ak+1]), I2(f ; [ak , ak+1])}
verwenden.
Numerische Mathematik I 192
Praktische Aspekte der Integration Typ 2: a-posteriori Fehlerschatzung
4.4.2 Typ 2: a-posteriori Fehlerschatzung
Es wird die Fehlerdarstellung der summierten Quadraturformel zur Schrittweiteh = (b − a)/N in der Form
∫ b
a
f (x) dx − IΣn (f ;N) = cf hk + r(f , h)hk+1
mit einer Konstanten cf und einer beschrankten Funktion r(f , .) vorausgesetzt.
Es stehen zwei Werte IΣn (f ;N1) und IΣn (f ;N2) zu den Schrittweiten h1 > h2 > 0,hj = (b − a)/Nj , j = 1, 2, zur Verfugung.
Dann kann die Konstante cf geschatzt werden durch Losen von
cf (hk1 − hk2 ) = −
(IΣn (f ;N1)− IΣn (f ;N2)
)+O(hk+1
1 ).
Fur h1 ≪ 1 folgt
cf ≈ − IΣn (f ;N1)− IΣn (f ;N2)
hk1 − hk2.
Die “richtige” Schrittweite h fur 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 tretenin den Teilintervallen ganz unterschiedliche Fehlergroß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 uberall halbiert, sondern nur in Intervallen[ak , ak+1], wo ein Fehlerschatzer signalisiert, dass
∣∣∣∣∫ ak+1
ak
f (x) dx − In(f , [ak , ak+1])
∣∣∣∣ >ǫ(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. DieIntegration ∫ 1
0
√|x − 0.7| dx
zur Genauigkeit ǫ = 10−4 und Ausgabe der einzelnen Unterteilungsintervalle (mit Hilfe dertrace-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 erhalt dieZwischenpunkte
(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 einewillkurliche 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:
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9adaptive summierte Simpson−Regel
Numerische Mathematik I 196