Numerická integrace. 6. listopadu 2012

Numerická integrace Michal Čihák

6. listopadu 2012

Výpočty integrálů v praxi

• V přednáškách z matematické analýzy jste se seznámili s mnoha

metodami výpočtu integrálů. • V praxi se ale poměrně často můžeme setkat s případy, kdy žádná

z těchto metod nevede k cíli. • Potom je jedinou možností použít některou z přibližných metod výpočtu integrálu.










Příklad: V matematické statistice často pracujeme s tzv. normovaným normálním rozdělením, jehož hustota je dána funkcí x2 1 f (x) = √ e− 2 2π



Chceme-li určit pravděpodobnost, že hodnota náhodné veličiny s normálním rozdělením leží v intervalu ha, bi, pak musíme vypočítat Z a

b

x2 1 √ e− 2 dx. 2π



Chceme-li určit pravděpodobnost, že hodnota náhodné veličiny s normálním rozdělením leží v intervalu ha, bi, pak musíme vypočítat Z a

b

x2 1 √ e− 2 dx. 2π

Tento integrál ale nelze exaktními metodami určit.

Dělené diference

Představme si, že známe hodnoty neznámé spojité funkce f v n + 1 bodech x0 < x1 < · · · < xn . Pro těchto n + 1 bodů existuje n + 1 tzv. nultých dělených diferencí funkce f f [xi ] = f (xi ),

i = 0, 1, . . . , n.

Dělené diference

Představme si, že známe hodnoty neznámé spojité funkce f v n + 1 bodech x0 < x1 < · · · < xn . Pro těchto n + 1 bodů existuje n + 1 tzv. nultých dělených diferencí funkce f f [xi ] = f (xi ),

i = 0, 1, . . . , n.

Dále existuje n tzv. prvních dělených diferencí funkce f f [xi , xi+1 ] =

f [xi+1 ] − f [xi ] , xi+1 − xi

i = 0, 1, . . . , n − 1.

Dělené diference Představme si, že známe hodnoty neznámé spojité funkce f v n + 1 bodech x0 < x1 < · · · < xn . Pro těchto n + 1 bodů existuje n + 1 tzv. nultých dělených diferencí funkce f f [xi ] = f (xi ),

i = 0, 1, . . . , n.

Dále existuje n tzv. prvních dělených diferencí funkce f f [xi , xi+1 ] =

f [xi+1 ] − f [xi ] , xi+1 − xi

i = 0, 1, . . . , n − 1.

Dále pokračujeme indukcí. Pokud známe (k − 1)-ní dělené diference f [xi , xi+1 , . . . , xi+k−1 ],

f [xi+1 , xi+2 , . . . , xi+k ],

potom pro k-tou dělenou diferenci platí f [xi , xi+1 , . . . , xi+k ] =

f [xi+1 , xi+2 , . . . , xi+k ] − f [xi , xi+1 , . . . , xi+k−1 ] . xi+k − xi

Dělené diference

Celý proces ukončíme určením jediné n-té dělené diference f [x0 , . . . , xn ] =

f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ] . xn − x0

Dělené diference

Celý proces ukončíme určením jediné n-té dělené diference f [x0 , . . . , xn ] =

f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ] . xn − x0

S pomocí dělených diferencí lze Lagrangeův interpolační polynom pro funkci f s uzly x0 , x1 , . . . , xn vyjádřit ve tvaru Pn (x) =f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 ) + . . . f [x0 , x1 , . . . , xn ](x − x0 )(x − x1 ) · · · (x − xn−1 ), který se nazývá Newtonův vzorec dělených diferencí.

Základní metody výpočtu určitých integrálů

• Funkci, jejíž určitý integrál na intervalu ha, bi chceme určit,

nahradíme Lagrangeovým interpolačním polynomem. • Z tohoto polynomu určíme určitý integrál na intervalu ha, bi. • Otázkou je, jaký stupeň Lagrangeova interpolačního polynomu zvolit

(kolik uzlů zvolit).









Obdélníkové pravidlo

Začneme tím, že zvolíme jeden uzel (stupeň Lagrangeova interpolačního polynomu bude 0). Tento uzel zvolíme uprostřed intervalu ha, bi.

Obdélníkové pravidlo Začneme tím, že zvolíme jeden uzel (stupeň Lagrangeova interpolačního polynomu bude 0). Tento uzel zvolíme uprostřed intervalu ha, bi.


Potom lze vyjádřit Z

b

Z f (x)dx ≈

a

b

Z P0 (x)dx =

a

b

f [x0 ]dx = f [x0 ](b−a) = f

a

a+b 2

(b−a).


Obdélníkové pravidlo: Z

b

f (x)dx ≈ f

a

a+b 2

(b − a).

Lichoběžníkové pravidlo

Nyní zvolíme dva uzly (stupeň Lagrangeova interpolačního polynomu bude 1). Za uzly zvolíme krajní body intervalu ha, bi.

Lichoběžníkové pravidlo Nyní zvolíme dva uzly (stupeň Lagrangeova interpolačního polynomu bude 1). Za uzly zvolíme krajní body intervalu ha, bi.



b

Z f (x)dx ≈

a

b

Z P1 (x)dx =

a

a

b

f [x0 ] + f [x0 , x1 ](x − x0 ) dx.


Postupně vypočítáme Z

b

Z f (x)dx ≈

a

b

Z P1 (x)dx =

a

b

f [x0 ] + f [x0 , x1 ](x − x0 ) dx =

a

b (x − a)2 = = f [a]x + f [a, b] 2 a f (b) − f (a) (b − a)2 (a − a)2 = f (a)(b − a) + − = b−a 2 2 f (a) + f (b) = (b − a) . 2


Lichoběžníkové pravidlo: Z

b

f (x)dx ≈ (b − a) a

f (a) + f (b) . 2

Simpsonovo pravidlo

Nyní zvolíme tři uzly (stupeň Lagrangeova interpolačního polynomu bude 2). Za uzly zvolíme krajní body a střed intervalu ha, bi.

Simpsonovo pravidlo Nyní zvolíme tři uzly (stupeň Lagrangeova interpolačního polynomu bude 2). Za uzly zvolíme krajní body a střed intervalu ha, bi.

Simpsonovo pravidlo


b

Z f (x)dx ≈

a

b

P2 (x)dx. a

Simpsonovo pravidlo

Simpsonovo pravidlo

Simpsonovo pravidlo: Z b (b − a) a+b f (x)dx ≈ f (a) + 4f + f (b) . 6 2 a

Porovnání jednotlivých metod

V první tabulce jsou uvedeny hodnoty určitých integrálů různých funkcí vypočtené na intervalu h1; 1,2i různými metodami – obdélníkovým pravidlem (Midpoint), lichoběžníkovým pravidlem (Trapezoidal) a Simpsonovým pravidlem. V prvním řádku tabulky jsou přitom uvedeny přesné hodnoty.

Porovnání jednotlivých metod V první tabulce jsou uvedeny hodnoty určitých integrálů různých funkcí vypočtené na intervalu h1; 1,2i různými metodami – obdélníkovým pravidlem (Midpoint), lichoběžníkovým pravidlem (Trapezoidal) a Simpsonovým pravidlem. V prvním řádku tabulky jsou přitom uvedeny přesné hodnoty.

Porovnání jednotlivých metod

Ve druhé tabulce jsou uvedeny hodnoty určitých integrálů různých funkcí vypočtené na intervalu h0; 2i různými metodami – obdélníkovým pravidlem (Midpoint), lichoběžníkovým pravidlem (Trapezoidal) a Simpsonovým pravidlem. V prvním řádku tabulky jsou přitom uvedeny přesné hodnoty.

Porovnání jednotlivých metod Ve druhé tabulce jsou uvedeny hodnoty určitých integrálů různých funkcí vypočtené na intervalu h0; 2i různými metodami – obdélníkovým pravidlem (Midpoint), lichoběžníkovým pravidlem (Trapezoidal) a Simpsonovým pravidlem. V prvním řádku tabulky jsou přitom uvedeny přesné hodnoty.

Jak zvýšit přesnost numerické integrace

Příklad: Určete pomocí Simpsonova pravidla

R2 0

ex dx.


Příklad: Určete pomocí Simpsonova pravidla Z 0

2

ex dx ≈

R2 0

ex dx.

1 0 (e + 4e1 + e2 ) = 6,4207278. 3


Příklad: Určete pomocí Simpsonova pravidla Z

2

ex dx ≈

0

R2 0

ex dx.

1 0 (e + 4e1 + e2 ) = 6,4207278. 3

Přesná hodnota je přitom Z 2 ex dx = [ex ]20 = e2 − e0 = 6,3890561. 0


Příklad: Určete pomocí Simpsonova pravidla Z

2

ex dx ≈

0

R2 0

ex dx.

1 0 (e + 4e1 + e2 ) = 6,4207278. 3

Přesná hodnota je přitom Z 2 ex dx = [ex ]20 = e2 − e0 = 6,3890561. 0

Absolutní chyba aproximace je tedy 0,0316717. Takováto chyba může být pro některé aplikace nepřijatelně vysoká.


Zkusme zvýšit přesnost aproximace tím, že rozdělíme interval h0, 2i na dva podintervaly h0, 1i a h1, 2i a na každém z nich použijeme Simpsonovo pravidlo


Zkusme zvýšit přesnost aproximace tím, že rozdělíme interval h0, 2i na dva podintervaly h0, 1i a h1, 2i a na každém z nich použijeme Simpsonovo pravidlo Z 2 Z 1 Z 2 1 1 ex dx = ex dx + ex dx ≈ (e0 + 4e0,5 + e1 ) + (e1 + 4e1,5 + e2 ) = 6 6 0 0 1 1 0 = (e + 4e0,5 + 2e1 + 4e1,5 + e2 ) = 6,3912102. 6


Zkusme zvýšit přesnost aproximace tím, že rozdělíme interval h0, 2i na dva podintervaly h0, 1i a h1, 2i a na každém z nich použijeme Simpsonovo pravidlo Z 2 Z 1 Z 2 1 1 ex dx ≈ (e0 + 4e0,5 + e1 ) + (e1 + 4e1,5 + e2 ) = ex dx + ex dx = 6 6 0 1 0 1 0 = (e + 4e0,5 + 2e1 + 4e1,5 + e2 ) = 6,3912102. 6 Absolutní chyba aproximace se zmenšila na 0,0021541, což je méně než 10 % původní absolutní chyby.


Zkusme ještě dále rozdělit intervaly h0, 1i a h1, 2i na další podintervaly. Při použití Simpsonova pravidla dostaneme

Jak zvýšit přesnost numerické integrace Zkusme ještě dále rozdělit intervaly h0, 1i a h1, 2i na další podintervaly. Při použití Simpsonova pravidla dostaneme Z 2 Z 0,5 Z 1 Z 1,5 Z 2 ex dx = ex dx + ex dx + ex dx + ex dx ≈ 0

0

0,5

1

1,5

1 1 0 (e + 4e0,25 + e0,5 ) + (e0,5 + 4e0,75 + e1 )+ ≈ 12 12 1 1 + (e1 + 4e1,25 + e1,5 ) + (e1,5 + 4e1,75 + e2 ) = 12 12 1 0 = (e + 4e0,25 + 2e0,5 + 4e0,75 + e1 + 4e1,25 + 2e1,5 + 4e1,75 + e2 ) 12 = 6,3891937.

Jak zvýšit přesnost numerické integrace Zkusme ještě dále rozdělit intervaly h0, 1i a h1, 2i na další podintervaly. Při použití Simpsonova pravidla dostaneme Z 2 Z 0,5 Z 1 Z 1,5 Z 2 x x x x e dx = e dx + e dx + e dx + ex dx ≈ 0

0

0,5

1

1,5

1 1 0 (e + 4e0,25 + e0,5 ) + (e0,5 + 4e0,75 + e1 )+ ≈ 12 12 1 1 + (e1 + 4e1,25 + e1,5 ) + (e1,5 + 4e1,75 + e2 ) = 12 12 1 0 = (e + 4e0,25 + 2e0,5 + 4e0,75 + e1 + 4e1,25 + 2e1,5 + 4e1,75 + e2 ) 12 = 6,3891937. Absolutní chyba aproximace se zmenšila na 0,0001376, což je už jen 0,4 % původní absolutní chyby (při použití Simpsonova pravidla na celý interval h0, 2i).

Složené Simpsonovo pravidlo

Zobecněním předchozího postupu získáme tzv. Simpsonovo složené pravidlo.

Složené Simpsonovo pravidlo Zobecněním předchozího postupu získáme tzv. Simpsonovo složené pravidlo.

Složené Simpsonovo pravidlo

Zvolíme sudé číslo n a rozdělíme interval ha, bi na n podintervalů. Označíme-li h = (b − a)/n, potom krajní body podintervalů jsou a = x0 < x1 < · · · < xn = b, kde xi = x0 + ih pro každé i = 0, 1, . . . , n. Na každém z intervalů (x2j−2 , x2j ), kde j = 1, 2, . . . , n/2, potom použijeme Simpsonovo pravidlo

Složené Simpsonovo pravidlo Zvolíme sudé číslo n a rozdělíme interval ha, bi na n podintervalů. Označíme-li h = (b − a)/n, potom krajní body podintervalů jsou a = x0 < x1 < · · · < xn = b, kde xi = x0 + ih pro každé i = 0, 1, . . . , n. Na každém z intervalů (x2j−2 , x2j ), kde j = 1, 2, . . . , n/2, potom použijeme Simpsonovo pravidlo Z

b

f (x)dx ≈ a

n/2 Z X j=1

=

x2j

f (x)dx =

x2j−2

n/2 X h j=1

3

[f (x2j−2 ) + 4f (x2j−1 ) + f (x2j )] =

 n/2 (n/2)−1 X X h = f (x0 ) + 2 f (x2j ) + 4 f (x2j−1 ) + f (xn ) . 3 j=1 j=1 

Absolutní chyba aproximace Simpsonovým složeným pravidlem

Lze odvodit (metodami diferenciálního a integrálního počtu), že absolutní chyba aproximace Simpsonovým složeným pravidlem je rovna 4 h (b − a) (4) , f (ξ) 180 kde f (4) (ξ) je čtvrtá derivace funkce f v bodě ξ, přičemž ξ je nějaké číslo z intervalu (a, b).

Složené Simpsonovo pravidlo – příklad

Rπ Příklad: Určete pomocí složeného Simpsonova pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat?


Rπ Příklad: Určete pomocí složeného Simpsonova pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 4 h (b − a) (4) h4 π h4 π π5 180 f (ξ) = 180 sin ξ ≤ 180 · 1 = 180n4 < 0,00002.


Rπ Příklad: Určete pomocí složeného Simpsonova pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 4 h (b − a) (4) h4 π h4 π π5 180 f (ξ) = 180 sin ξ ≤ 180 · 1 = 180n4 < 0,00002. Z poslední nerovnosti určíme, že n > 18. Můžeme tedy zvolit například n = 20 a h = π/20.

Složené Simpsonovo pravidlo – příklad Rπ Příklad: Určete pomocí složeného Simpsonova pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 4 h (b − a) (4) h4 π h4 π π5 = f (ξ) sin ξ ≤ ·1= < 0,00002. 180 180 180 180n4 Z poslední nerovnosti určíme, že n > 18. Můžeme tedy zvolit například n = 20 a h = π/20. S použitím těchto hodnot obdržíme pomocí složeného Simpsonova pravidla   Z π 9 10 X X π  jπ (2j − 1)π sin 0 + 2 +4 + sin π  = sin xdx ≈ sin sin 60 10 20 0 j=1 j=1 = 2,00000678.

Složené Simpsonovo pravidlo – příklad Rπ Příklad: Určete pomocí složeného Simpsonova pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 4 h (b − a) (4) h4 π h4 π π5 180 f (ξ) = 180 sin ξ ≤ 180 · 1 = 180n4 < 0,00002. Z poslední nerovnosti určíme, že n > 18. Můžeme tedy zvolit například n = 20 a h = π/20. S použitím těchto hodnot obdržíme pomocí složeného Simpsonova pravidla   Z π 9 10 X X π  jπ (2j − 1)π sin 0 + 2 sin +4 sin + sin π  = sin xdx ≈ 60 10 20 0 j=1 j=1 = 2,00000678. Rπ Přesná hodnota je přitom 0 sin xdx = 2, absolutní chyba je tedy v tomto případě rovna 0,00000678, což je skutečně méně než zadaná maximální přípustná chyba 0,00002.

Složené Simpsonovo pravidlo – shrnutí

Předpokládejme, že funkce f má spojité derivace až do 4. řádu na intervalu ha, bi. Nechť dále n je sudé číslo, h = (b − a)/n a xj = a + jh pro každé j = 0, 1, . . . , n. Potom pro nějaké ξ ∈ (a, b) platí   Z b (n/2)−1 n/2 X X h f (x)dx = f (a) + 2 f (x2j ) + 4 f (x2j−1 ) + f (b) − 3 a j=1 j=1 −

h4 (b − a) (4) f (ξ). 180

Složené lichoběžníkové pravidlo

Předpokládejme, že funkce f má spojité derivace až do 2. řádu na intervalu ha, bi. Nechť dále n je libovolné přirozené číslo, h = (b − a)/n a xj = a + jh pro každé j = 0, 1, . . . , n. Potom pro nějaké ξ ∈ (a, b) platí   Z b (n−1 X h h2 (b − a) (2) f (x)dx = f (a) + 2 f (xj ) + f (b) − f (ξ). 2 12 a j=1

Složené obdélníkové pravidlo

Předpokládejme, že funkce f má spojité derivace až do 2. řádu na intervalu ha, bi. Nechť dále n je sudé číslo, h = (b − a)/(n + 2) a xj = a + (j + 1)h pro každé j = −1, 0, 1, . . . , n + 1. Potom pro nějaké ξ ∈ (a, b) platí Z

b

f (x)dx = 2h a

n/2 X j=0

f (x2j ) +

h2 (b − a) (2) f (ξ). 6

Složené lichoběžníkové pravidlo – příklad

Rπ Příklad: Určete pomocí složeného lichoběžníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat?


Rπ Příklad: Určete pomocí složeného lichoběžníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 2 h (b − a) (2) h2 π h2 π π3 f (ξ) = sin ξ ≤ ·1= < 0,00002. 12 12 12 12n2


Rπ Příklad: Určete pomocí složeného lichoběžníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 2 h (b − a) (2) h2 π h2 π π3 = f (ξ) sin ξ ≤ ·1= < 0,00002. 12 12 12 12n2 Z poslední nerovnosti určíme, že n > 359 (oproti Simpsonovu pravidlu výrazně více!)

Složené lichoběžníkové pravidlo – příklad Rπ Příklad: Určete pomocí složeného lichoběžníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 2 h (b − a) (2) h2 π h2 π π3 f (ξ) = sin ξ ≤ ·1= < 0,00002. 12 12 12 12n2 Z poslední nerovnosti určíme, že n > 359 (oproti Simpsonovu pravidlu výrazně více!) Pro srovnání ještě uveďme, že pro n = 20 a h = π/20 obdržíme pomocí složeného lichoběžníkového pravidla Z π sin xdx ≈ 1,995886. 0

Absolutní chyba této aproximace je 0,004114, což jen potvrzuje naše předchozí zjištění, že 20 intervalů je pro dosažení požadované přesnosti málo.

Složené obdélníkové pravidlo – příklad

Rπ Příklad: Určete pomocí složeného obdélníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat?


Rπ Příklad: Určete pomocí složeného obdélníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 2 h (b − a) (2) h2 π h2 π π3 f (ξ) = sin ξ ≤ ·1= < 0,00002. 6 6 6 6(n + 2)2


Rπ Příklad: Určete pomocí složeného obdélníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 2 h (b − a) (2) h2 π h2 π π3 = f (ξ) sin ξ ≤ ·1= < 0,00002. 6 6 6 6(n + 2)2 Z poslední nerovnosti určíme, že n + 2 > 508 (oproti Simpsonovu pravidlu opět výrazně více!)

Složené obdélníkové pravidlo – příklad Rπ Příklad: Určete pomocí složeného obdélníkového pravidla 0 sin dx s absolutní chybou menší než 0,00002. Kolik podintervalů intervalu h0, πi budeme pro tento účel potřebovat? 2 h (b − a) (2) h2 π h2 π π3 f (ξ) = sin ξ ≤ ·1= < 0,00002. 6 6 6 6(n + 2)2 Z poslední nerovnosti určíme, že n + 2 > 508 (oproti Simpsonovu pravidlu opět výrazně více!) Pro srovnání ještě uveďme, že pro n + 2 = 20 a h = π/20 obdržíme pomocí složeného obdélníkového pravidla Z π sin xdx ≈ 2,008248. 0

Absolutní chyba této aproximace je 0,008248, což jen potvrzuje naše předchozí zjištění, že 20 intervalů je pro dosažení požadované přesnosti málo.

Numerická integrace. 6. listopadu 2012

Recommend Documents