Kap4 Numeriche Mathe
Kap4 Numeriche Mathe
Kap4 Numeriche Mathe
Numerische Integration
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
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
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.
b−a
ak = a + k · , k = 0, 1, . . . , N,
N
ak−1 +ak
und verwenden der Mittelpunkte yk = 2 , also
Systematik zur Definition der “einfachen” Regeln: Gegeben sei das Integral
Rb
a f (x) dx.
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
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
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
Die Exaktheit für Polynome vom Grad ≤ n führt zur Darstellung des Fehlers der
Quadraturformel.
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.
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
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
Die Simpsonregel mit N = 4 Teilintervallen liefert schon 4 Nachkommastellen. Sie erfordert die
Auswertung des Integranden an 9 Stellen.
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
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 .
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 .
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
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]
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
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
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
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
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
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
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)
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
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
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.
Die Positivität der Gewichte der Gauß-Formeln bewirkt, dass die Operator-Norm
von In auf C [−1, 1]
X N−1
h h
a(h) = I1Σ (f ; N) = f (a) + h f (a + jh) + f (b).
2 2
j=1
c) Rekursionsformel:
X
k−1
k Bj
B0 = 1, Bk = − für k = 1, 2, . . .
j=0
j k −j +1
(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)
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)!
Das feste Vorzeichen von b2m+2 (x) − B2m+2 ergibt sich aus Eigenschaft (ii).
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].
Für 0 ≤ k ≤ m gilt
Z b
2k+2
f (x) dx − aj,k = O(hj−k ) für k ≤ j → ∞.
a
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
gilt.
verwenden.
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 .
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)
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]
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