Metody Numeryczne, semestr letni 2013/2014
Wykład 9: Całkowanie numeryczne, cz.2
Kwadratury Gaussa.
b
n
Przypomnijmy: całkę oznaczoną S( f ) = R p( x) f ( x) dx przybliżamy sumą Q( f ) = P Akf ( xk), a
k=0
xk ∈ [ a, b].
Za kryterium dokładności kwadratury przyjmujemy zgodność Q( W ) i S( W ), gdy W jest wielomianem.
Definicja 1 Mówimy, że kwadratura Q jest rzędu r, r 1 , gdy 1. S( W ) = Q( W ) dla wszystkich wielomianów W ( x) stopnia mniejszego niż r, 2. istnieje wielomian W ( x) stopnia r taki, że S( W ) 6= Q( W ) .
Rozpatrzmy problem wyboru węzłów i współczynników przy ustalonej wadze p( x) i liczbie n tak, aby rząd kwadratury był jak najwyższy. Taką kwadraturę nazywamy kwadraturą Gaussa.
n
Twierdzenie 1 Rząd kwadratury Q( f ) = P Akf ( xk) wynosi co najmniej n + 1 wtedy i tylko k=0
wtedy, gdy kwadratura ta jest postaci
n
Q( f ) = S( L
X
n) =
Akf ( xk) ,
k=0
gdzie
b
Z
( x − x
A
0)( x − x 1) . . . ( x − xk− 1)( x − xk+1) . . . ( x − xN ) k =
p( x)Φ k( x) dx,
Φ k( x) =
.
( xk − x 0)( xk − x 1) . . . ( xk − xk a
− 1)( xk − xk+1) . . . ( xk − xN ) Na podstawie powyższego twierdzenia wnioskujemy, że problem kwadratury Gaussa spro-wadza się do odpowiedniego wyboru węzłów xk, k = 0 , 1 , . . . , n.
Kwadratury Gaussa są związane z wielomianami ortogonalnymi.
Twierdzenie 2 Wielomiany ortogonalne mają tylko pierwiastki rzeczywiste jednokrotne leżące w przedziale ( a, b) .
n
Twierdzenie 3 Nie istnieje kwadratura postaci Q( f ) = P Akf ( xk) rzędu wyższego niż 2( n + 1) .
k=0
n
b
Twierdzenie 4 Kwadratura Q( f ) = P A R
k f ( xk ) o współczynnikach Ak =
p( x)Φ k( x) dx jest k=0
a
rzędu 2( n + 1) wtedy i tylko wtedy, gdy xk są pierwiastkami wielomianu ortogonalnego Pn+1( x) , czyli wielomianu ortogonalnego stopnia n + 1 z wagą p w przedziale ( a, b) .
Twierdzenie 5 Wszystkie współczynniki Ak w kwadraturach Gaussa są dodatnie.
Zauważmy, że kwadratura Gaussa jest rzędu 2 n + 2, czyli jest dokładna dla wielomianów stopnia 2 n + 1. Kwadratury Gaussa zależą od wagi p.
1
Kwadratury Gaussa dla przedziału skończonego 1) Kwadraturę najwyższego rzędu z wagą p( x) ≡ 1 nazywamy kwadraturą Gaussa-Legendre’a.
Przyjmujemy [ a, b] = [ − 1 , 1]. Dla wagi p( x) = 1 ciąg wielomianów ortogonalnych tworzą wielomiany Legendre’a:
1
dn
Pn( x) =
( x 2 − 1) n.
2 nn! dxn
Współczynniki i oszacowanie błędu kwadratur Gaussa-Legendre’a są następujące: 2
Ak = −
,
k = 0 , 1 , . . . , n
( n + 2) Pn+2( xk) P ′n+1( xk) 22 n+3(( n + 1)!)4
|E( f) | ¬
M
(2 n + 3)((2 n + 2)!)3
2 n+2 ,
gdzie x
k są miejscami zerowymi wielomianu Legendre’a Pn+1( x) oraz M 2 n+2 =
max f (2 n+2)( x) .
x∈[ − 1 , 1]
Tablica początkowych węzłów i współczynników kwadratur Gaussa-Legendre’a na przedziale [ − 1 , 1]: n k
xk
Ak
1
0
− 0 , 577350
1
1
0,577350
1
2
0
− 0 , 774597 0,555556
1
0
0,888889
2
0,774597
0,555556
3
0
− 0 , 861136 0,347855
1
− 0 , 339981 0,652145
2
0,339981
0,652145
3
0,861136
0,347855
4
0
− 0 , 906180 0,236927
1
− 0 , 538469 0,478629
2
0
0,568889
3
0,538469
0,478629
4
0,906180
0,236927
Uwaga 1 Dla dowolnego przedziału całkowania [ a, b] w całce dokonujemy liniowej zamiany zmiennych:
b + a
b − a
t =
+
x,
x ∈ [ − 1 , 1] ,
2
2
czyli
b
1
Z
b − a Z
f ( t) dt =
g( x) dx,
2
a
− 1
gdzie g( x) = f b+ a + b−ax .
2
2
Do tak przekształconej całki stosujemy podane wzory i dostajemy b
n
n
Z
b − a
b − a
f ( t) dt ≈ Q( f) =
X A
X A
2
k f ( tk ) =
2
kg( xk ) ,
a
k=0
k=0
2
gdzie xk - miejsca zerowe wielomianu Legendre’a Pn+1( x) oraz b + a
b − a
tk =
+
x
2
2
k
Uwaga 2 Wagę p( x) ≡ 1 wybieramy wtedy, gdy funkcja podcałkowa jest gładka na przedziale
[ a, b].
Uwaga 3 Oszacowanie błędu wzoru Gaussa-Legendre’a na przedziale [ a, b] wynosi ( b − a)2 n+3(( n + 1)!)4
|E( f) | ¬
M
f (2 n+2)( x) .
(2 n + 3)((2 n + 2)!)3
2 n+2 ,
M 2 n+2 = max
x∈[ a,b]
Przykład 1 Stosując wzór Gaussa-Legendre’a oparty na czterech węzłach ( n = 3) obliczyć 1 √
całkę R
1 + 2 t 3 dt.
0
Mamy a = 0, b = 1, więc dokonujemy liniowej zamiany zmiennych, która sprowadzi do całkowania po [ − 1 , 1]:
t = 1 + 1 x
1
2
2
s
3
Z
√
3
dt = 1 dx
1 Z 1
1
1
1
1 + 2 t 3 dt =
X
2
=
1 + 2
+ x
dx ≈
A
t
0
1
2
2
2
2
k g( xk) ,
− 1
0
k=0
x
− 1 1
gdzie
s
s
1
1
3
1
g( xk) =
1 + 2
+ x
=
1 +
(1 + x
2
2 k
4
k )3 ,
k = 0 , 1 , 2 , 3 .
Obliczamy
q
k
xk
g( xk) =
1 + 1 (1 + x
4
k )3
Ak
Akg( xk)
0
-0,861136
1,000335
0,347855
0,347971
1
-0,339981
1,035316
0,652145
0,675176
2
0,339981
1,265504
0,652145
0,825292
3
0,861136
1,616064
0,347855
0,562156
3
2,410596= P Akg( xk)
k=0
Ostatecznie
1
1 3
Z
√
X A
1 + 2 t 3 dt
2
k g( xk) = 1 , 205298 ≈
k=0
0
Oszacowanie błędu wynosi
(1 − 0)9(4!)4
(4!)4 · 94
|E( f) | ¬
max f (8)( t) ¬
≈ 5 , 3 · 10 − 8 .
9 · (8!)3
t∈[0 , 1]
9 · (8!)3
Funkcja podcałkowa może posiadać osobliwości polegające na tym, że jest ona nieograniczo-na na [ a, b] lub jej pochodne są nieograniczone. Taką funkcję przedstawiamy w postaci iloczynu p · f, gdzie p zawiera wszelkie osobliwości funkcji podcałkowej, a f jest już funkcją dostatecznie gładką na przedziale [ a, b]. Przyjęcie funkcji p za funkcję wagową umożliwia usunięcie części osobliwej zarówno z sumy Q( f ), jak i z błędu E( f ). Oczywiście zakładamy, że p jest nieujemną, całkowalną funkcją na [ a, b].
3
1
√
na przedziale [ − 1 , 1] nazywamy kwadraturami Gaussa-1 −x 2
Czebyszewa. Ciąg wielomianów ortogonalnych tworzą wielomiany Czebyszewa (pierwszego rodzaju):
Tn( x) = cos ( n arc cos x) .
Współczynniki, węzły i oszacowanie błędu są następujące: π
Ak =
,
k = 0 , 1 , . . . , n
n + 1
(2 k + 1) π
xk = cos 2 n + 2 π
|E( f) | ¬
M
f (2 n+2)( x) .
22 n+1(2 n + 2)!
2 n+2 ,
M 2 n+2 = max
x∈[ − 1 , 1]
2
Przykład 2 Obliczyć całkę R
t 2
√
dt stosując kwadraturę Gaussa-Czebyszewa dla n = 2.
0
t(2 −t)
Wprowadzamy nową zmienną, która należy do przedziału [ − 1 , 1]: t = 1 + x
2
1
2
Z
t 2
dx = dt
Z
(1 + x)2
dt =
X
=
√
dx ≈
A
q
k f ( xk) =
t(2
t
0
2
1 − x 2
0
− t)
k=0
− 1
x
− 1 1
π
π 2
3 π 2
5 π 2!
=
1 + cos
+ 1 + cos
+ 1 + cos
3
6
6
6
√
√
π
3!2
3!2
=
1 +
+ (1 + 0)2 + 1 −
= 1 , 5 π ≈ 4 , 712389
3
2
2
gdzie
f ( xk) = (1 + xk)2 ,
k = 0 , 1 , 2 .
4
Kwadratury Gaussa dla przedziału nieskończonego 3) Na przedziale ( −∞, + ∞) przyjmujemy wagę p( x) = e−x 2, dla której ciąg wielomianów ortogonalnych tworzą wielomiany Hermite’a:
Hn( x) = ( − 1) nex 2 dn e−x 2 .
dxn
Kwadraturę nazywamy w tym przypadku kwadraturą Gaussa-Hermite’a.
Współczynniki, węzły i oszacowanie błędu:
2 n+2( n + 1)!
Ak =
,
k = 0 , 1 , . . . , n
Hn+2( xk) H′n+1( xk)
√
( n + 1)! π
|E( f) | ¬
M
f (2 n+2)( x) .
2 n+1(2 n + 2)!
2 n+2 ,
M 2 n+2 = max
x
∈ R
gdzie xk są miejscami zerowymi wielomianu Hermite’a Hn+1( x). Tablica początkowych węzłów i współczynników kwadratur Gaussa-Hermite’a:
n k
xk
Ak
1
0
− 0 , 707107 0,886227
1
0,707107
0,886227
2
0
− 1 , 224745 0,295409
1
0
1,181636
2
1,224745
0,295409
3
0
− 1 , 650680 0,081313
1
− 0 , 524648 0,804914
2
0,524648
0,804914
3
1,650680
0,081313
4
0
− 2 , 020183 0,019953
1
− 0 , 958572 0,393619
2
0
0,945309
3
0,958572
0,393619
4
2,020183
0,019953
4) Dla wagi p( x) = e−x i przedziału [0 , + ∞) ciąg wielomianów ortogonalnych tworzą wielomiany Laguerre’a:
dn
Ln( x) = ( − 1) nex
xne−x .
dxn
Kwadraturę nazywamy kwadraturą Gaussa-Laguerre’a.
Współczynniki, węzły i oszacowanie błędu:
(( n + 1)!)2
Ak =
,
k = 0 , 1 , . . . , n
Ln+2( xk) L′n+1( xk)
(( n + 1)!)2
|E( f) | ¬
M
f (2 n+2)( x) .
(2 n + 2)!
2 n+2 ,
M 2 n+2 = max
x∈[0 , + ∞
gdzie xk są miejscami zerowymi wielomianu Laguerre’a Ln+1( x).
Tablica początkowych węzłów i współczynników kwadratur Gaussa-Laguerre’a: 5
xk
Ak
1
0
0,585786
0,853553
1
3,414214
0,146447
2
0
0,415775
0,711093
1
2,294280
0,278518
2
6,289945
0,010389
3
0
0,322548
0,603154
1
1,745761
0,357419
2
4,536620
0,038888
3
9,395071
0,000539
4
0
0,263560
0,521756
1
1,413403
0,398667
2
3,596426
0,075942
3
7,085810
0,003612
4
12,640801 0,000032
+ ∞
Przykład 3 Obliczyć całkę R e−xx 3 dx stosując kwadraturę Gaussa-Laguerre’a dla n = 2.
0
+ ∞
Z
e−xx 3 dx ≈ 0 , 711093 ·(0 , 415775)3+0 , 278518 ·(2 , 294280)3+0 , 010389 ·(6 , 289945)3 ≈ 5 , 999938 .
0
Oszacowanie błędu wynosi E( f ) = 0, ponieważ dla funkcji f ( x) = x 3 jest max |f 6( x) | = 0.
x∈ R
Dla porównania wartość dokładna całki wynosi 6.
6