1.1 1.1.1
Numerické integrování Úvodní úvahy
Naším cílem bude přibližný numerický výpočet určitého integrálu Zb
I=
f (x)dx.
(1.1)
a
Je-li známa k integrované funkci f primitivní funkce F (F 0 (x) = f (x)), můžeme integrál v (1.1) spočítat analyticky, Zb
f (x)dx = F (b) − F (a).
(1.2)
a
Tedy např.: ·
Zπ
¸π
= 2.
sin(x)dx = − cos x
(1.3)
0
0
Obvykle ale primitivní funkci neznáme a integrál (1.1) musíme počítat numericky. Mluvíme o numerické kvadratuře. Lze např. využít součtové definice určitého integrálu. Interval ha, bi rozdělíme body x0 ≡ a < x1 < x2 < . . . < xn−1 < xn ≡ b
(1.4)
na dostatečně malé intervaly hxk−1 , xk i, k = 1, 2, . . . , n. Pak I≈
n X
f (ξk )∆xk ,
(1.5)
k=1
kde ∆xk = xk − xk−1 je šířka k-tého podintervalu a ξk je jeho libovolný bod, např. ξk = xk nebo ξk = xk−12+xk . Pro „rozumnéÿ funkce f konverguje suma na pravé straně vztahu (1.5) pro ∆xk → 0 k přesné hodnotě I. Jiná metoda spočívá v nahrazení funkce f (x) vhodnou aproximující funkcí p(x), jejíž integrál dokážeme spočítat analyticky. Z přibližné rovnosti f (x) ≈ p(x) vyplývá Zb
I≈
p(x)dx.
(1.6)
a
Funkci p(x) volíme zpravidla ve formě interpolačního polynomu. Polynom n-tého stupně p(x) = qn xn + qn−1 xn−1 + . . . + q1 x + q0
(1.7)
obsahuje n + 1 koeficientů qn , qn−1 , . . . , q1 , q0 . Ty jsou určeny soustavou n + 1 rovnic p(xk ) = f (xk ), k = 0, 1, 2, . . . , n.
Po rozepsání xn0 qn + xn−1 qn−1 + . . . + x0 q1 + q0 = f (x0 ), 0 xn1 qn + xn−1 qn−1 + . . . + x1 q1 + q0 = f (x1 ), 1 .. .
(1.8)
xnn qn + xn−1 n qn−1 + . . . + xn q1 + q0 = f (xn ). Dělící body xk intervalu ha, bi (viz(1.2)) jsou zpravidla voleny jako ekvidistantní. V praxi se obvykle používá kombinace obou předchozích metod - součtové a interpolační: interval ha, bi se rozdělí na podintervaly hxk−1 , xk i, Zb
I=
f (x)dx =
N X
Zxk
f (x)dx,
(1.9)
k=1 xk−1
a
a poté se funkce f (x) aproximuje funkcí p(x) zvlášť na jednotlivých podintervalech, Zxk
Zxk
f (x)dx ≈ xk−1
(1.10)
p(x)dx. xk−1
Jsou-li podintervaly hxk−1 , xk i dostatečně úzké, vystačíme s polynomy p nízkého stupně.
1.1.2
Interpolační metoda
Naším úkolem je na intervalu ha, bi interpolovat funkci f (x) polynomem n-tého stupně p(x) (1.7). Interpolační uzly xk , k = 0, 1, 2, . . . , n zvolíme pro jednoduchost jako ekvidistantní, b−a . (1.11) n Omezíme se na polynom nultého stupně (aproximace funkce f (x) konstantou, obdélníková metoda interpolace), prvního stupně (aproximace lineární funkcí, lichoběžníková metoda) a stupně druhého (aproximace kvadratickou funkcí, Simpsonova metoda). xk = a + kh, k = 0, 1, 2, . . . , n, h ≡
a) Obdélníková metoda Obdélníková metoda vychází z aproximace f (x) ≈ q0 pro x ∈ ha, bi,
(1.12)
kde q0 je konstanta. Zvolíme-li q0 = f ( a+b ), pak 2 Zb
I=
f (x)dx ≈ a
Ã
Zb
f a
!
a+b dx = (b − a)f 2
Ã
!
a+b . 2
(1.13)
Geometricky odpovídá tato aproximace nahrazení plochy pod křivkou f (x) plochou obdélníka o výšce f ( a+b ) (obr.1). 2
b) Lichoběžníková metoda V tomto případě je funkce f (x) na intervalu ha, bi aproximována funkcí lineární, f (x) ≈ p(x) ≡ q1 x + q0 .
(1.14)
Koeficienty p1 , p0 určíme z interpolačních podmínek p(a) = f (a), p(b) = f (b): aq1 + q0 = f (a), bq1 + q0 = f (b).
(1.15)
Řešením této jednoduché soustavy je f (b) − f (a) bf (a) − af (b) , q0 = . b−a b−a Aproximace integrálu je tedy q1 =
Zb
I≈ a
"
x2 (q1 x + q0 )dx = q1 + q0 x 2
#b
= (b − a) a
f (a) + f (b) . 2
(1.16)
(1.17)
Geometricky odpovídá tato aproximace nahrazení plochy pod křivkou f (x) lichoběžníkem o výšce b − a a základnách f (a), f (b) (obr.2).
c) Simpsonova metoda Funkce f je na intervalu ha, bi aproximována kvadratickou funkcí p(x) = q2 x2 + q1 x + q0 . Interval ha, bi rozdělíme uzly a+b , x2 = b (1.18) 2 na dva podintervaly o šířce h = b−a . Koeficienty q2 , q1 , q0 jsou určeny soustavou tří 2 lineárních rovnic (1.8) pro n = 2. Tento postup je dosti těžkopádný, a to pracujeme s polynomy pouze 2. řádu! Elegantnější způsob nalezení interpolačního polynomu nabízí Lagrangeova metoda. Polynom p(x) n-tého řádu, procházející n + 1 body [xk , f (xk )], k = 0, 1, 2, . . . , n, je v této metodě vyjádřen ve tvaru x0 = a, x1 =
p(x) =
n X
f (xi )li (x),
(1.19)
i=1
kde li (x) =
(x − x0 )(x − x1 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn ) (xi − x0 )(xi − x1 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )
(1.20)
jsou tzv. Lagrangeovy polynomy. Ty jsou zkonstruovány tak, aby *
li (xk ) =
1, 0,
k=i k 6= i,
(1.21)
takže p(xk ) =
n X
(1.22)
f (xi )li (xk ) = f (xk )1,
i=0
tj. graf polynomu p(x) prochází body [xk , f (xk )]. Z přibližného vztahu f (x) ≈ p(x) dostáváme Zb
I≈
p(x)dx =
Zb X n
f (xi )li (x)dx =
a i=0
a
n X
Zb
f (xi )wi , wi ≡
i=0
li (x)dx.
(1.23)
a
Koeficienty wi se nazývají váhy v uzlech xi . Speciálně pro polynom 2. řádu s uzly (1.18), dává Lagrangeova metoda interpolaci f (x) ≈ p(x) = f (x0 ) +f (x1 )
(x − x1 )(x − x2 ) + (x0 − x1 )(x0 − x2 )
(x − x0 )(x − x2 ) (x − x0 )(x − x1 ) + f (x2 ) . (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )
(1.24)
Při výpočtu váhových faktorů w0 , w1 , w2 budeme využívat substituci x = a + ht, dx = hdt, 0 ≤ t ≤ 2. Zb
w0 = a
2
Z (x − x1 )(x − x2 ) (a + ht − a − h)(a + ht − a − 2h) dx = hdt = (x0 − x1 )(x0 − x2 ) (−h)(−2h) 0
2
1 hZ = (t − 1)(t − 2)dt = h, 2 3 0
{z
|
}
2 3
Zb
w1 = a
2
Z (x − x0 )(x − x2 ) (a + ht − a)(a + ht − a − 2h) dx = hdt = (x1 − x0 )(x1 − x2 ) h(−h) 0
Z2
= −h 0
Zb
w2 = a
=
h 2
4 t(t − 2)dt = h, 3 2
Z (a + ht − a)(a + ht − a − h) (x − x0 )(x − x1 ) dx = hdt = (x2 − x0 )(x2 − x1 ) 2hh 0
Z2 0
1 t(t − 1)dt = h. 3
(1.25)
Simpsonův vzorec má tedy tvar 1 4 1 h I ≈ f (x0 ) h + f (x1 ) h + f (x2 ) h = [f (x0 ) + 4f (x1 ) + f (x2 )] 3 3 3 3 Analogicky bychom odvodili interpolační formule vyšších řádů.
1.1.3
(1.26)
Složené vzorce
Čím užší bude interval integrace, tím přesnější bude aproximace funkce f polynomem p(x). Interval ha, bi proto rozložíme na malé podintervaly hxk−1 , xk i, na nichž lze chybu aproximace f (x) ≈ p(x) očekávat relativně malou. Uzly xk ze vztahu (1.4) zvolíme pro jednoduchost ekvidistantní s krokem h, xk = a + kh, k = 0, 1, 2, . . . , n, h =
b−a . n
(1.27)
a) Složená obdélníková metoda Formule (1.9) spolu s formulí (1.13) aplikovanou na dílčí intervaly hxk−1 , xk i dávají pro tento případ · µ
I≈h f
¶
µ
¶
µ
x0 + x1 x1 + x2 xn−1 + xn +f + ... + f 2 2 2
¶¸
(1.28)
.
b) Složená lichoběžníková metoda Lichoběžníkovou aproximaci (1.17) aplikujeme na jednotlivé intervaly hxk−1 , xk i: h I≈ 2
·µ
¶
µ
¶
µ
¶¸
f (x0 ) + f (x1 ) + f (x1 ) + f (x2 ) + . . . + f (xn−1 ) + f (xn )
,
(1.29)
tj. ·
¸
1 1 I ≈ h f (x0 ) + f (x1 ) + f (x2 ) + . . . + f (xn−1 ) + f (xn ) . 2 2
(1.30)
c) Složená Simpsonova metoda V tomto případě musíme mít sudý počet uzlů, n = 2m, h = b−a . Simpsonův vzorec 2m (1.26) aplikujeme postupně na intervaly hx0 , x2 i, hx2 , x4 i, . . . , hx2m−1 , x2m i: h I≈ 3
·µ
¶
µ
¶
f (x0 ) + 4f (x1 ) + f (x2 ) + f (x2 ) + 4f (x3 ) + f (x4 ) + µ
¶¸
+ . . . + f (x2m−2 ) + 4f (x2m−1 ) + f (x2m )
.
(1.31)
Po úpravě ·
h I≈ f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 )+ 3 ¸ + . . . + 2f (x2m−2 ) + 4f (x2m−1 ) + f (x2m ) .
(1.32)
Povšimněte si, že i složené kvadraturní vzorce lze vyjádřit pomocí tabelizovaných hodnot f (xk ) a váhových faktorů wk v uzlech xk , Zb
f (x)dx ≈
n X
f (xk )wk .
(1.33)
k=0
a
Pro lichoběžníkovou metodu je w0 = wn =
h , w1 = w2 = . . . = wn−1 = h, 2
(1.34)
kdežto pro Simpsonovu metodu 1 4 w0 = w2m = h, w1 = w3 = . . . = w2m−1 = h, 3 3 2 w2 = w4 = . . . = w2m−2 = h. 3
1.1.4
(1.35)
Cykly složených metod, Richardsonova metoda
Označme numerickou kvadraturu integrálu symbolem N a chybu integrace jako E. Přesná hodnota I je tedy I = N + E.
(1.36)
Pro lichoběžníkovou metodu lze chybu E vyjádřit ve tvaru (b − a) 00 f (η)h2 , (1.37) 12 kde f 00 (η) značí druhou derivaci integrované funkce f v blíže nespecifikovaném bodě η intervalu ha, bi. Lichoběžníková metoda poskytuje tedy přesný výsledek pro lineární funkce, které mají druhou derivaci nulovou. E=−
Protože bod η není znám, odhadujeme chybu E shora. Označme Mk = max |f (k) (x)|. ha,bi
Pak (b − a) M2 h2 . 12 Chyba je druhého řádu v mocninách kroku h. Lze ji rovněž vyjádřit ve formě |E| ≤
E = Ch2 + členy vyšších řádů v h,
(1.38)
(1.39)
kde C je konstanta nezávislá na kroku h. Pro složenou Simpsonovu metodu lze odvodit (b − a) (1.39) f (η)h4 . (1.40) 180 Simpsonova metoda dává přesnou hodnotu integrálu pro polynomy do třetího řádu, neboť mají čtvrtou derivaci nulovou. E=−
Analogií vztahů (1.39), (1.40) jsou pro Simpsonovu metodu vztahy |E| ≤
(b − a) M4 h4 , 180
E = Ch4 + členy vyšších řádů.
(1.41) (1.42)
Požadované přesnosti integrace lze dosáhnout zmenšováním integračního kroku. Mámeli navíc numerickou integraci provedenou pro dva různé kroky h1 , h2 , můžeme provést tzv. Richardsonovu extrapolaci na krok h = 0, odpovídající přesné hodnotě integrálu. Pro kroky h1 , h2 pišme . I = N1 + Chk1 , . I = N2 + Chk2 ,
(1.43)
V tomto vztahu jsou N1 , N2 numerické hodnoty integrálu získané s kroky h1 , h2 . V rozvoji chyby E (vztahy (1.39), (1.42)) jsme se omezili na hlavní člen, přičemž k = 2 pro lichoběžníkovou a k = 4 pro Simpsonovu metodu. Vztahy (1.43) jsou proto jen přibližné. Vztahy (1.43) představují soustavu dvou rovnic pro neznámé I, C. Aproximaci přesné hodnoty I dostaneme nejsnáze tak, že první rovnici vynásobíme hk2 , druhou hk1 , a obě rovnice poté od sebe odečteme. Dostaneme k k . N2 h 1 − N1 h 2 I= . hk1 − hk2
(1.44)
Při zjemňování kroku h je výhodné v každé iteraci zdvojnásobit počet dělících bodů, neboť pak v následujícím kroku plně využijeme dělící body a funkční hodnoty z kroku předchozího. Dosadíme-li h1 = 2h2 do (1.44), obdržíme k . 2 N2 − N1 . I= 2k − 1 Speciálně pro lichoběžníkovou metodu dává předchozí formule
(1.45)
. 4N2 − N1 I= , (1.46) 3 to však odpovídá metodě Simpsonově. Rozdělme interval ha, bi na n = 2m stejných dílků s rozestupem h a uzly xk ,
b−a , xk = a + kh, k = 0, 1, 2, . . . , 2m. 2m Kvadratura s krokem h dává h=
µ
¶
1 1 N2 = h y0 + y1 + y2 + . . . + y2M , 2 2 kdežto kvadratura s dvojnásobným krokem dává µ
(1.48)
¶
1 1 N1 = 2h y0 + y2 + y4 + . . . + y2M . 2 2 Po dosazení (1.48), (1.49) do (1.46) dostaneme . h I = (y0 + 4y1 + 2y2 + . . . + y2M ) , 3 což je ale Simpsonova sumace.
1.1.5
(1.47)
(1.49)
(1.50)
Jak je to v Matlabu
Funkce QUAD QU AD = Numerické vyhodnocení integrálu, adaptivní Simpsonova metoda. Q = QU AD(F U N, A, B) aproximuje integral funkce F U N v mezích od A do B, kromě chyby 1.e − 6 s použitím rekurzivní adaptivní Sipsonovy metody. Funkce Y = F U N (X) pracuje s vektorem X a jako výsledek vrací vektor Y , vyhodnocený integrand každého z prvků vektoru X. Q = QU AD(F U N, A, B, T OL) používá absolutní chybu tolerance T OL namísto implicitní hodnoty 1.e − 6. Větší hodnoty tolerance T OL, poté proběhne méně výpočtů funkce a tedy rychleji výpočet, ale na úkor přesnosti výsledků. Funkce QU AD ve verzi MATLAB 5.3 používá méně spolehlivý algoritmus a implicitní tolerance má hodnotu 1.e − 3. [Q, F CN T ] = QU AD(. . .) vrací počet vyhodnocení funkce. QU AD(F U N, A, B, T OL, T RACE) s nenulovou stopou (TRACE) ukáže hodnotu [f cnt a b − a Q] během rekurze. QU AD(F U N, A, B, T OL, T RACE, P 1, P 2, ...) poskytuje jako další argumenty P 1, P 2, . . . předány přímo funkci F U N, F U N (X, P 1, P 2, . . .). Průchod prázdné matice co se týče T OL nebo T RACE k použití implicitních hodnot. Použití maticových operátorů .? , ./ a .∧ v definici F U N tak, že to lze vyhodnotit s vektorovým argumentem (vektor). Funkce QU ADL může pracovat více účinně s vysokou přesností a hladkou integrovanou funkcí.
Příklad: F U N můžeme zadat třemi různými způsoby. Řetězcové vyjádření umocnění jednoduché proměnné: Q = quad(0 1./(x.∧ 3 − 2? x − 5)0 , 0, 2); Inline objekt: F = inline(0 1./(x.∧ 3 − 2? x − 5)0 ); Q = quad(F, 0, 2); Ukazatel funkce (operátor @): Q = quad(@myfun, 0, 2); kde myfun.m je M-soubor: function y = myfun(x) y = 1./(x.∧ 3 − 2? x − 5);
Funkce QUADL QU ADL = číselné vyhodnocení integrálu, adaptivní Lobattova metoda. Q = QU ADL(F U N, A, B) vyšetřujeme přibližnou hodnotu integrálu funkce F U N v mezích od A do B, kromě chyby 1.e − 6 s použitím vysoce uspořádané rekurzivní adaptivní metody. Funkce Y = F U N (X) pracuje s vektorem X a jako výsledek vrací vektor Y , vyhodnocený integrand každého z prvků vektoru X. Q = QU ADL(F U N, A, B, T OL) používá absolutní chybu tolerance T OL namísto implicitní hodnoty 1.e − 6. Větší hodnoty tolerance T OL, poté proběhne méně výpočtů funkce a tedy rychleji výpočet, ale na úkor přesnosti výsledků. [Q, F CN T ] = QU ADL(. . .) vrací počet vyhodnocení funkce. QU ADL(F U N, A, B, T OL, T RACE) s nenulovou stopou (TRACE) ukáže hodnotu [f cnt a b − a Q] během rekurze. QU ADL(F U N, A, B, T OL, T RACE, P 1, P 2, . . .) poskytuje jako další argumenty P 1, P 2, . . . předány přímo funkci F U N , F U N (X, P 1, P 2, . . .). Průchod prázdné matice co se týče T OL nebo T RACE k použití implicitních hodnot. Použití maticových operátorů .? , ./ a .∧ v definici F U N takže to lze vyhodnotit s vektorovým argumentem (vektor). Příklad: F U N můžeme zadat třemi různými způsoby.
Řetězcové vyjádření umocnění jednoduché proměnné: Q = quadl(0 1./(x.∧ 3 − 2? x − 5)0 , 0, 2); Inline objekt: F = inline(0 1./(x.∧ 3 − 2? x − 5)0 ); Q = quadl(F, 0, 2); Ukazatel funkce (operátor @): Q = quadl(@myfun, 0, 2); kde myfun.m je M-soubor: function y = myfun(x) y = 1./(x.∧ 3 − 2? x − 5);
