1
Numerische Integration durch Extrapolation
Pablo Thiel
1
Romberg-Verfahren
Idee: Im Gegensatz zur numerischen Integration mit Hilfe der einfachen bzw. zusammengesetzten Trapez-, Simpson-, 3/8- oder zum Beispiel der Milne-Regel, wird bei der Integration durch Extrapolation eine Ann¨aherung an den Integralwert durch eine Extrapolation u ¨ber diskrete Werte, die eine zusammengesetzte Regel mit verschiedenen Schrittweiten liefert, erreicht. Das heisst, dass das gesuchte Integral I mit Hilfe einer zusammengesetzen Regel f¨ ur ¨ verschiedene Schrittweiten h approximiert wird. Uber diese gefundenen Werte wird nun ein Interpolationspolynom gelegt, und an der Stelle h = 0 ausgewertet. Da 0 ausserhalb des Intervalls liegt, in dem das Interpolationspolynom entstand, nennt man das Verfahren Extrapolation. Man extrapoliert bei h = 0, da anzunehmen ist, dass wenn man immer kleiner werdende Schrittweiten betrachtet, die Approximation immer genauer wird. Die Annahme also, dass lim = I ist. h→0
Zur numerischen Integration des Integrals
Rb
f (x) dx durch Extrapolation verwenden
a
wir die Trapezregel T (h) zusammen mit der Romberg-Folge hi 2i −1
X hi (f (a + jhj )) + f (b)) T (hi ) = (f (a) + 2 2 j=0 mit der Romberg-Folge hi :=
b−a 2i
f¨ ur i = 1, 2, ..., m.
2
Numerische Integration durch Extrapolation
Dabei kann man T (hi+1 ) rekursiv aus T (hi ) berechnen. Es gilt: 1 1 T (hi+1 ) = T ( hi ) = T (hi ) + hi+1 (f (a + hi+1 ) + f (a + 3hi+1 ) + ... + f (b − hi+1 )) 2 2
T (h) besitzt folgende Entwicklung welche in Abschnitt 2 bewiesen wird. T (h) = τ0 + τ1 h2 + τ2 h4 + ... + τm h2m + αm+1 (h)h2m+2 mit τ0 :=
Rb
f (x) dx. Dabei sind die τi von h unabh¨angige Konstanten und αm+1 (h) ist
a
eine beschr¨ankte Funktion.
Da man T (h) nur f¨ ur h > 0 nicht aber f¨ ur h = 0 berechnen kann, interpoliert man T an den St¨ utzstellen h0 , h1 , ..., hm > 0 mit den Funktionswerten T (h0 ), T (h1 ), ..., T (hm ) und extrapoliert f¨ ur h = 0 mit Hilfe des Neville-Schemas. Neville-Schema zur Polynomwertbestimmung: xi x0 x1 x2 . . . . xm
Ti,0 T0,0 T1,0 T2,0
Tm,0
= T (xi ) T1+i,1 = T (x0 ) T1,1 = T (x1 ) T2,1 = T (x2 ) T3,1 . . . . . . . Tm,1 = T (xm )
T2+i,2 T2,2 T3,2 T4,2 . . Tm,2
... Tm−1+i,m−1 ... Tm−1,m−1 ... Tm,m−1 ...
Tm,m Tm,m
mit m ≥ i ≥ k ≥ 1 und Ti,0 = T (hi ) Ti,k (x) =
(x−xi−k )Ti,k−1 −(x−xi )Ti−1,k−1 (x) xi −xi−k
k≥1
Daraus folgt: (x−xi )Ti,k−1 +(xi −xi−k )Ti,k−1 −(x−xi )Ti−1,k−1 (x) xi −xi−k T −Ti−1,k−1 Ti,k−1 + i,k−1 x−xi−k −1 x−x
Ti,k = =
i
¨ (Der Ubersicht halber schreibe ich Ti,k statt Ti,k (x), da x fest ist.)
1 Romberg-Verfahren
3
Da hier an der Stelle x = 0 ausgewertet werden soll und T (h) eine Entwicklung in h2 hat muss also xi = h2i gesetzt werden. Es gilt also der sogenannte Neville-Aitken-Algorithmus:
Ti,k (0) = Ti,k−1 +
Ti,k−1 − Ti−1,k−1 hi−k 2 hi
−1
Hieraus kann man verschiedene interessante Zusammenh¨ange ableiten: Zum Beispiel erh¨alt man f¨ ur T1,1 die Simsonregel: T1,0 −T0,0 4 T1,1 = T1,0 + h0 2 = 3 T1,0 − 13 T0,0 = =
−1 h1 b−a (f (a) + 2f ( a+b ) + f (b)) − b−a (f (a) + f (b)) 3 2 6 1 a+b (b − a)(f (a) + 4f ( ) + f (b)) (Simpsonregel) 6 2
Analog folgt zum Beispiel auch bei T2,2 die Milne-Regel und f¨ ur hi = pulcherima“zur¨ uck. ”
b−a 3i
in T1,1 die
Obrige Rechnung l¨asst sich auch generell f¨ ur Ti,1 machen: Ti,1 = 34 Ti,0 − 13 Ti−1,0 = 43 T (hi ) − 13 T (2hi ) = 34 (τ0 + τ1 h2i + τ2 h4i + ... + τm h2m + ...) i − 13 (τ0 + τ1 4h2i + τ2 16h4i + ... + τm 22m h2m + ...) i 4 6 = τ0 − 4τ2 hi − 20τ3 hi ... Das heisst also, dass der Fehler in der zweiten Spalte des Neville-Schemas wie h4i gegen 0 strebt. Insgesammt strebt der Fehler in der k-ten Spalte (d.h. von Ti,k−1 ) wie h2k i gegen 0.
Beispiel: Gesucht ist eine N¨aherung f¨ ur
R2 1
xi 1 1 2 1 4 1 8 1 16
1 x
dx = log(2) = 0.6931471806...
Ti,0 = T (xi ) T1+i,1 T2+i,2 T3+i,3 T4,4 ¯ 0.75 0.694 0.6931746032 0.6931474776 0.6931471819 0.78¯3 0.6932539683 0.6931479015 0.69314791831 0.6970238095 0.6931545307 0.6931471843 0.6941218504 0.6931476528 0.6933912022
Unterstrichen sind hier jeweils die mit dem exaktem Ergebniss u ¨bereinstimmenden Stellen.
4
2
Numerische Integration durch Extrapolation
Euler-Maclaurinsche Summenformel
Satz 2.1: F¨ ur f ∈ ζ 2m+2 [a, b] (d.h. f ist auf [a, b] (2m + 2) mal stetig differenzierbar) besitzt die Trapezsumme die Entwicklung T (h) = τ0 + τ1 h2 + τ2 h4 + ... + τm h2m + αm+1 (h)h2m+2 mit τ0 :=
Rb
(1)
f (x) dx. Dabei sind die τi von h unabh¨angige Konstanten und αm+1 (h) ist
a
eine beschr¨ankte Funktion.
Beweis: Zum Beweis ben¨otigt man zuerst folgende Vorbemerkungen: Bernulli Polynome: Bk (x) mit k = 0, 1, ... rekursiv definiert durch: (a) B0 (x) = 1 0 (b) Bk (x) = kBk−1 (x) k ≥ 1 R1 (c) Bk (x) dx = 0 k≥1
(2)
0
Aus (2) folgt: (d) Bk (x) = Ak +
Rx
Bk−1 (t) dt k ≥ 1
0
mit der Konstanten Ak = Bk (0) so, dass
R1
Bk (x) dx = 0 .
0
Insbesondere ist z.B.: B1 (x) = x − 12 B2 (x) = x2 − x + 16 B3 (x) = x3 − 32 x2 + 12 x B4 (x) = x4 − 2x3 + x2 −
1 30
F¨ ur Bk (x) gelten die folgenden wichtigen Eigenschafften: (a) Bk (0) = Bk (1) (b) Die Polynome Pk (x) := Bk ( 21 + x) sind f¨ ur gerades k gerade, f¨ ur ungerades k ungerade Polynome in x. (c) B2k+1 (0) = B2k+1 (1) = 0
k≥2 (3) ∀k ≥ 1
2 Euler-Maclaurinsche Summenformel
5
Beweis: zu a) Aus (2b) und (2c) folgt: Bk (0) − Bk (1) = Ak + k
R0
Bk−1 (t) dt − Ak − k
0
R1
Bk−1 (t) dt = 0
0
zu b) P0 (x) = B0 ( 21 + x) = 1 ist gerade. Sei nun P2k (x) f¨ ur ein k ≥ 0 ein gerades Polynom. Es gilt weiterhin wegen (2b) und (2c): 0 0 P2k+1 (x) = B2k+1 (x + 12 ) = (2k + 1)B2k (x + 12 ) = (2k + 1)P2 k(x) 1 R2 R1 R1 und P2k+1 (x) dx = P2k+1 (x − 12 ) dx = B2k+1 (x) dx = 0 0
− 12
0
Da P2k (x) gerade ist, ist das Polynom Rx q(x) := (2k + 1) P2k (t) dt ein ungerades Polynom und es gilt: 0 − 12
1 2
R
q(x) dx = −
− 12
R
− 12
q(−x) dx =
1 2
R
1
q(x) dx = −
1 2
R2
q(x) dx = 0
− 12 1
0
Damit hat q(x) die charakteristischen Eigenschaften (q (x) = (2k+1)P2k (x);
R2
q(x) dx =
− 12
0) von P2k+1 (x): ⇒ q(x) = P2k+1 (x) ⇒ P2k+1 (x) ist ungerade. Rx Wegen P2k+2 = (2k + 2) P2k+1 (t) dt + A mit A konstant ist P2k+2 (x) eine gerade Funk0
tion. zu c) Aus (3a) und (3b) folgt: B2k+1 (1) = B2k+1 (0) = P2k+1 (− 21 ) = −P2k+1 ( 12 ) = −B2k+1 (1) = 0.
Als Bernullische Zahlen bezeichnent man: Bk := (−1)k+1 B2k (0) Also zum Beispiel: B1 = 61 , B2 =
1 30
, B3 =
1 42
, B4 =
1 . 30
6
Numerische Integration durch Extrapolation
Nun definieren wir uns 1-periodische Funktionen: Sk : R → R mit: Sk (x) := Bk (x − i) mit x ∈ I = {x|i ≤ x < i + 1} , i ganz.
Diese gen¨ ugen folgender Rekursion: Sk (x) = Sk (0) + k
Rx
Sk+1 (t) dt k ≥ 2
(4)
0
da wegen (2d) gilt: R R Sk (x) = Bk (x − i) = Ak + k −0x−i Bk−1 (t) dt = Bk (0) + k −0x−i Sk−1 (t + i) dt x−1 R = Sk (0) + k Sk−1 (t) dt f¨r k ≥ 2 0
Ausserdem gilt wegen (3) f¨ ur alle k ≥ 1: (a) Bk (0) = Sk (0) = Sk (1) = Sk (2) = ... (b) S2k+1 (0) = 0 (c) S2k (0) = (−1)k+1 Bk .
Beweis: zu a) gilt, da f¨ ur alle x ∈ Z gilt: Sk (x) = Sk (0) + k
Rx
Sk−1 (t) dt = Sk (0) = Bk (0)
0
zu b) S2k+1 (0) = B2k+1 (0) = 0 zu c) S2k (0) = B2k (0) = (−1)k+1 Bk
Zum Beweis von Satz 2.1 betrachten wir statt f die Funktion g mit g(t) := f (a + th)
(5)
2 Euler-Maclaurinsche Summenformel
7
Durch schrittweiser partieller Integration von
Rn
0
S1 (t)g (t) dt folgt:
0
Zn
g(0) g(n) S1 (t)g (t) dt = + g(1) + ... + g(n − 1) + − 2 2 0
g(t) dt
(6)
Bk (2k−1) (g (n) − g (2k−1) (0)) + Rm+1 (2k)!
(7)
0
Man kann
Zn 0
Rn
0
S1 (t)g (t) dt aber auch so umformen:
0
Zn
0
S1 (t)g (t) dt =
m X
(−1)k+1
k=1
0
mit dem Restglied Rm+1
1 = (2m + 2)!
Zn
(S2m+2 (t) − S2m+2 (0)) = g (2m+2) (t) dt
(8)
0
Wegen g(t) = f (a + th) folgt durch einfaches umrechnen: Rn
(a)
0 (k)
(b) (c)
g(t) dt =
g
1 h
Rb
a k (k)
(t) = h f
f (x) dx (a + th)
k = 0, 1, 2, ...
(9)
g(0) g(n) +g(1)+...+g(n−1)+ 2 = 2 f (a) f (b) 1 +f (a+h)+...+f (b−h)+ 2 = h T (h) 2
¨ Aus diesen Uberlegungen und nach (6),(7) und (8) gilt nun: T (h) = τ0 + τ1 h2 + τ2 h4 + ... + τm h2m + αm+1 (h)h2m+2 wegen: T (h) = h2 (f (a) + f (a + h) + ... + f (b − h) + f (b)) = h( g(0) + g(1) + ... + g(n − 1) + g(n) ) 2 2 n m R P Bk (g (2k+1) (n) − g (2k−1) (0))) + Rm+1 ) = h( g(t) dt + ((−1)k+1 (2k)! 0
=
Rb
=
Rb
f (x) dx + h
k=1
a
a
k=1 m P
f (x) dx +
Pm
k+1
( (−1)(2k)! Bk (h2k−1 f (2k−1) (a + tn) − h2k−1 f (2k−1) (a + 0n))) + hRm+1
k=1 (
(−1)k+1 Bk (f (2k−1) (b) (2k)!
− f (2k−1) (a))) + hRm+1
8
Numerische Integration durch Extrapolation
Rn −1 Mit dem Restglied: hRm+1 = h (2m+2)! (S2m+2 (t) − S− 2m + 2(0))g (2m+2) (t) dt 0
Rb −1 = h (2m+2)! h2m+1 (S2m+2 ( x−a ) − S2m+2 (0))f (2m+2) (x) dx h a
−1 = h2m+2 (2m+2)!
Rb a
f (2m+2) (x)(S2m+2 ( s−a ) − S2m+2 (0)) dx h
Genau erh¨alt man also die Euler-Maclaurinsche Summenformel durch: τ0 =
Rb
f (x)dx
a
τk = αm+1 (h) =
(−1)k+1 Bk (f (2k−1) (b) (2k)!
−1 (2m+2)!
Rb a
− f (2k−1) (a)
k = 1, 2, ..., m
f (2m+2) (x)(S2m+2 ( x−a ) − S2m+2 (0)) dx h
S2m+2 ist eine stetige periodische Funktion, es gibt also eine von h unabh¨angige obere Schranke f¨ ur |αm+1 (h)|. Damit ist Satz 2.1 bewiesen.
