integrály lze vypočítat snadno pomocí tabulek a klasických integračních metod jako je per partes nebo substituce. Tak například integrály

Numerické metody

6.

´ INTEGROVAN ´ Í A DERIVOVAN ´ Í 6. NUMERICKE

Numerick´ e integrov´ an´ı a derivov´ an´ı

Pr˚ uvodce studiem Pˇri ˇreˇsen´ı r˚ uzných u ´ loh je potˇreba nalézt hodnotu urˇcitého integrálu. Nˇekteré integrály lze vypoˇc´ıtat snadno pomoc´ı tabulek a klasických integraˇcn´ıch metod jako je per partes nebo substituce. Tak napˇr´ıklad integrály π 1 x e dx a cos x dx 0

0

ˇ maj´ı hodnotu e − 1 a 0. Casto vˇsak staˇc´ı malá zmˇena v integrované funkci a klasický v´ ypoˇcet se stane sloˇzit´ ym nebo dokonce nemoˇzným. Pokuste se vypoˇc´ıtat integrály

1

x2

e dx

a

0

π

cos x2 dx

0

a uvid´ıte jak daleko se vám podaˇr´ı doj´ıt! Klasické integraˇcn´ı metody nelze pouˇz´ıt v˚ ubec, jestliˇze integrovaná funkce je dána tabulkou (napˇr. to m˚ uˇze být v´ ysledek mˇeˇren´ı nebo pˇredchoz´ıho v´ ypoˇctu). V této kapitole se budeme zabývat metodami numerick´ ymi pro obecnˇe zadanou u ´ lohu: pˇredpokládáme, ˇze je dána spojitá funkce f na intervalu a, b a máme vypoˇc´ıtat hodnotu urˇcitého integrálu b f (x) dx. I=

(6.0.1)

a

Z geometrického pohledu pˇredstavuje ˇc´ıslo I velikost plochy obrazce, kter´ y je vymezen grafem funkce f . Numerické metody jsou navrhovány tak, ˇze poˇc´ıtaj´ı velikost plochy pˇribliˇzného obrazce pomoc´ı nˇekolika funkˇcn´ıch hodnot. Lze je tedy pouˇz´ıt pro jakoukoliv funkci f , dokonce i pro funkci danou tabulkou, a jejich pracnost je relativnˇe malá. Pˇri numerickém v´ ypoˇctu derivace je situace podobná.

- 113 -


6.1.

Numerické metody

Newton-Cotesovy vzorce

C´ıle Odvod´ıme Newton-Cotesovy vzorce a ukáˇzeme jejich pouˇzit´ı.

Pˇredpokl´ adan´ e znalosti Interpolaˇcn´ı polynom. Interpolaˇcn´ı chyba.

V´ yklad Naˇs´ım c´ılem je navrhnout vzorce pro pˇribliˇzný v´ ypoˇcet integrálu (6.0.1). Ze základn´ıho kurzu integráln´ıho poˇctu v´ıme, ˇze snadno lze integrovat polynomy. Budeme proto postupovat tak, ˇze k funkci f sestav´ıme interpolaˇcn´ı polynom pn a ten integrujeme m´ısto f , tj. b b I= f (x) dx ≈ pn (x) dx. a

(6.1.1)

a

Pro jednoduchost uvaˇzujme na intervalu a, b ekvidistantn´ı uzly xi = a + iτ,

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

kde τ = (b−a)/n. Interpolaˇcn´ı polynom k funkci f m˚ uˇzeme zapsat v Lagrangeovˇe tvaru (viz odstavec 5.1.1), tj. pn (x) =

n

f (xi )ϕi (x),

kde ϕi (x) =

i=0

n $ x − xj . xi − xj j=0 j =i

Dosazen´ım v´ yrazu pro pn do pravé strany pˇribliˇzné rovnosti (6.1.1) dostaneme Newton-Cotesovy vzorce:

b

f (x) dx ≈ a

kde

n

wi f (xi ),

(6.1.2)

i=0

b

wi =

ϕi (x) dx,

i = 0, 1, . . . , n

a

- 114 -

(6.1.3)


Numerické metody

jsou integraˇcn´ı v´ ahy. Abychom lépe ilustrovali tento postup, odvod´ıme základn´ı integraˇcn´ı pravidla. a) Obd´ eln´ıkov´ e pravidlo. V tomto pˇr´ıpadˇe poloˇz´ıme n = 0, x0 =

a+b 2

a za

interpolaˇcn´ı polynom vezmeme konstantn´ı funkci p0 (x) = f ( a+b ). Proto 2

b a+b f (x) dx ≈ (b − a)f (6.1.4) = IObd . 2 a b) Lichobˇ eˇ zn´ıkov´ e pravidlo dostaneme pro n = 1 a x0 = a, x1 = b. Interpolaˇcn´ı polynom je pˇr´ımka p1 (x) = Po jej´ı integraci je

x−a x−b f (a) + f (b). a−b b−a

b

f (x) dx ≈ a

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

c) Simpsonovo pravidlo vznikne pro n = 2 a x0 = a, x1 =

(6.1.5) a+b , 2

x2 = b, kdy

je interpolaˇcn´ı polynom kvadratická funkce. Vypoˇctˇeme integraˇcn´ı váhy w0 , w1 a w2 . Dostáváme

b

w0 =

ϕ0 (x) dx =

a

a 1

= −1

s pouˇzit´ım substituce x = 1 (b 6

b

(x − x1 )(x − x2 ) dx (x0 − x1 )(x0 − x2 )

t(t − 1) b − a b−a dt = 2 2 6

b−a t + a+b . 2 2

Podobnˇe vypoˇc´ıtáme w1 = 46 (b − a) a w2 =

− a), takˇze

b a+b b−a [f (a) + 4f f (x) dx ≈ + f (b)] = ISimps . 6 2 a

Pˇ r´ıklad 6.1.1.

(6.1.6)

Pomoc´ı Newton-Cotesov´ ych vzorc˚ u (6.1.4), (6.1.5) a (6.1.6)

vypoˇctˇetˇe pˇribliˇznou hodnotu integrálu 1 ex dx. I= −1

- 115 -


Numerické metody

ˇ sen´ı: Máme a = −1, b = 1 a f (x) = ex . Podle obdéln´ıkového pravidla (6.1.4) Reˇ je IObd = 2e0 = 2. Lichobˇeˇzn´ıkové pravidlo (6.1.5) dává 2 ILich = [e−1 + e1 ] = 3.086161 2 a koneˇcnˇe pomoc´ı Simpsonova pravidla (6.1.6) vypoˇcteme 2 ISimps = [e−1 + 4e0 + e1 ] = 2.362054. 6 Poznamenejme jeˇstˇe, ˇze pˇresná hodnota integrálu je I = e1 − e−1 = 2.350402. Pozn´ amka Nejpˇresnˇejˇs´ı hodnoty jsme dos´ ahli pomoc´ı Simpsonova pravidla. To n´ as m˚ uˇze vést k domnˇence, ˇze jeˇstˇe pˇresnˇejˇs´ıch hodnot dos´ ahneme, pouˇzijeme-li NewtonCotesovy vzorce odvozené z interpolaˇcn´ıch polynom˚ u vyˇsˇs´ıch stupn˚ u. Obecnˇe to vˇsak nen´ı pravda. V Rungeho pˇr´ıkladu (pˇr´ıklad 5.1.4.) jsme vidˇeli,ˇze interpolaˇcn´ı polynom vysokého stupnˇe m˚ uˇze být velice ˇspatnou aproximac´ı. Newton-Cotesovy vzorce d´ avaj´ı v takovém pˇr´ıpadˇe nesmyslné výsledky. Následuj´ıc´ı vˇeta ukazuje vyjádˇren´ı integraˇcn´ı chyby. Vˇ eta 6.1.1. (i) Necht’ funkce f má na intervalu a, b spojitou druhou derivaci. Pak plat´ı: I − IObd I − ILich

f (ξ) (b − a)3 , = 24 f (ξ) (b − a)3 . = − 12

(6.1.7)

(ii) Necht’ funkce f má na intervalu a, b spojitou ˇctvrtou derivaci. Pak plat´ı: I − ISimps = −

f (4) (ξ) (b − a)3 . 90

Ve vˇsech pˇr´ıpadech je ξ bl´ıˇze neurˇcený bod z intervalu (a, b).

- 116 -

Numerické metody


D˚ ukaz: D˚ ukaz provedeme pouze pro lichobˇeˇzn´ıkové pravidlo. Vyjádˇren´ı interpolaˇcn´ı chyby z vˇety 5.1.2. má pro interpolaˇcn´ı polynom p1 s uzly x0 = a a x1 = b tvar:

˜ f (ξ) (x − a)(x − b). 2

f (x) − p1 (x) =

Integrac´ı (a uˇzit´ım vˇety o stˇredn´ı hodnotˇe integráln´ıho poˇctu) dostaneme f (ξ) b I − ILich = (x − a)(x − b) dx 2 a f (ξ) b−a f (ξ) = t(t + a − b) dx = − (b − a)3 . 2 12 0 Pouˇzili jsme substituci t = x − a.

2

Kontroln´ı ot´ azky Otázka 1. Jak vzniknou Newton-Cotesovy vzorce? Otázka 2. Jaké jsou základn´ı pravidla pro numerick´ y v´ ypoˇcet integrálu? Otázka 3. Znázornˇete graficky smysl obdéln´ıkového a lichobˇeˇzn´ıkového pravidla. Otázka 4. Proˇc nen´ı rozumné pouˇz´ıvat Newton-Cotesovy vzorce odvozené z interpolaˇcn´ıch polynom˚ u vysokého stupnˇe?

´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Pomoc´ı obdéln´ıkového, lichobˇeˇzn´ıkového a Simpsonova pravidla vypoˇctˇete pˇribliˇzné hodnoty integrál˚ u 1 2 a) ex dx;

π

b)

0

cos x2 dx.

0

V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. a) IObd = 1.284025, ILich = 1.859141, ISimps = 1.475731; b) IObd = 1.961189, ILich = −0.675916, ISimps = 1.082154.

- 117 -


6.2.

Numerické metody

Sloˇ zen´ e vzorce

C´ıle V pˇredchoz´ım odstavci jsme odvodili vzorce pro pˇribliˇzný v´ ypoˇcet interálu integrac´ı interpolaˇcn´ıho polynomu. Zároveˇ n jsme upozornili na to, ˇze vzorce odvozené z polynom˚ u vysokého stupnˇe mohou dávat nesmyslné v´ ysledky. Pro pˇresnˇejˇs´ı v´ ypoˇcty integrál˚ u se proto pouˇz´ıvaj´ı sloˇzené vzorce, které dostaneme slepen´ım” ” základn´ıch integraˇcn´ı pravidel.

Pˇredpokl´ adan´ e znalosti Obdéln´ıkové, lichobˇeˇzn´ıkové a Simpsonovo pravidlo. Chyby tˇechto pravidel.

V´ yklad Interval a, b rozdˇel´ıme na m stejnˇe dlouh´ ych d´ılk˚ u s krokem h = (b − a)/m, m ≥ 2, tj. pouˇzijeme uzly xi = a + ih = a +

i (b − a), m

i = 0, 1, . . . m.

Integrál rozloˇz´ıme následovnˇe:

b

I=

f (x) dx = a

m i=1

xi

f (x) dx.

(6.2.1)

xi−1

Jestliˇze kaˇzdý d´ılˇc´ı integrál nahrad´ıme pomoc´ı obdéln´ıkového nebo lichobˇeˇzn´ıkového pravidla, dostaneme odpov´ıdaj´ıc´ı pravidlo sloˇzené. a) Sloˇ zen´ e obd´ eln´ıkov´ e pravidlo jsme jiˇz odvodili v odstavci 1.1. Má tvar

! x1 + x2 xm−1 + xm x0 + x1 +f + ...+ f = h f 2 2 2

m xi−1 + xi = h . (6.2.2) f 2 i=1

ISO

- 118 -


Numerické metody

b) Sloˇ zen´ e lichobˇ eˇ zn´ıkov´ e pravidlo. V (6.2.1) provedeme náhradu podle vzorce (6.1.5), tj.

xi

f (x) dx ≈

xi−1

h [f (xi−1 ) + f (xi )], 2

a dostaneme ISL

1 1 = h f (x0 ) + f (x1 ) + . . . + f (xm−1 ) + f (xm ) 2 2 % & m−1 1 h f (x0 ) + 2 f (xi ) + f (xm ) . = 2 i=1

!

(6.2.3)

c) Sloˇ zen´ e Simpsonovo pravidlo. V tomto pˇr´ıpadˇe rozdˇel´ıme interval a, b na 2m stejnˇe dlouh´ ych d´ılk˚ u, m ≥ 2, s krokem h = (b − a)/(2m), takˇze budeme m´ıt lich´ y poˇcet uzl˚ u xi = a + ih, i = 0, 1, . . . 2m. Simpsonovo pravidlo (6.1.6) pouˇzijeme na kaˇzdém intervalu x2i−2 , x2i , i = 1, . . . , m, coˇz dává b m x2i m 2h I= f (x) dx = f (x) dx ≈ [f (x2i−2 ) + 4f (x2i−1 ) + f (x2i )]. 6 a i=1 x2i−2 i=1 Výraz na pravé stranˇe pˇrep´ıˇseme do pˇrehlednˇejˇs´ıho tvaru: h [f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + 2f (x4 ) + . . . + f (x2m )] 3 % & m m−1 h = f (x0 ) + 4 f (x2i−1 ) + 2 f (x2i ) + f (x2m ) . (6.2.4) 3 i=1 i=1

ISS =

Pˇ r´ıklad 6.2.1. Pomoc´ı pravidel (6.2.2), (6.2.3) a (6.2.4) vypoˇctˇetˇe pˇribliˇznou hodnotu integrálu

1

I=

ex dx.

−1

Interval −1, 1 pˇritom rozdˇelte na ˇctyˇri ˇcásti. ˇ sen´ı: Reˇ

(a) Poloˇz´ıme m = 4, takˇze krok bude h = 0.5 a dostaneme uzly

x0 = −1, x1 = −0.5, x2 = 0, x3 = 0.5 a x4 = 1. Sloˇzené obdéln´ıkového pravidlo (6.2.2) má tvar . ISO = 0.5[e−0.75 + e−0.25 + e0.25 + e0.75 ] = 2.326096.

- 119 -


Numerické metody

(b) Sloˇzené lichobˇeˇzn´ıkového pravidlo (6.2.3) vypoˇc´ıtáme pro stejné uzly, tj. ISL =

1 . · 0.5[e−1 + 2e−0.5 + 2e0 + 2e0.5 + e1 ] = 2.399166. 2

(c) U sloˇzeného Simpsonova pravidla (6.2.4) vezmeme m = 2, coˇz znamená, ˇze pouˇzijeme stejné dˇelen´ı intevalu jako v pˇredchoz´ıch pˇr´ıpadech (protoˇze h = 2/(2m) = 0.5). Dostaneme ISS =

0.5 −1 . [e + 4e−0.5 + 2e0 + 4e0.5 + e1 ] = 2.351195. 3

V následuj´ıc´ı vˇetˇe popisujeme ˇra´d jednotliv´ ych integraˇcn´ıch pravidel. Vˇ eta 6.2.1. (i) Necht’ má funkce f na intervalu a, b spojitou druhou derivaci. Sloˇzené obdéln´ıkové a sloˇzené lichobˇeˇzn´ıkové pravidlo jsou druhého ˇra´du, tj. plat´ı: |I − ISO | ≤ C1 h2 , |I − ISL | ≤ C2 h2 . (ii) Necht’ funkce f má na intervalu a, b spojitou ˇctvrtou derivaci. Sloˇzené Simpsonovo pravidlo je ˇctvrtého ˇra´du, tj. plat´ı: |I − ISS | ≤ C3 h4 . Konstanty C1 , C2 a C3 jsou bl´ıˇze neurˇcená nezáporná ˇc´ısla, která nezávis´ı na kroku h. D˚ ukaz: Omez´ıme se na lichobˇeˇzn´ıkové pravidlo. Staˇc´ı si uvˇedomit,ˇze jednoduché lichobˇeˇzn´ıkové pravidlo je ve sloˇzeném lichobˇeˇzn´ıkovém pravidle obsaˇzeno mkrát. Jestliˇze tedy chceme vyjádˇrit chybu, seˇcteme m-krát v´ yraz pro chybu jednoduchého pravidla (6.1.7). Dostaneme |I − ISL | = |

m i=1

f (ξi) f (ξi) 3 M 2 M(b − a) 2 h |≤ h3 | ≤ h h, − | h= 12 12 12 i=1 12 i=1 m

m

kde jsme pouˇzili M = maxζ∈a,b |f (ζ)|. Vid´ıme, ˇze C2 = M(b − a)/12.

- 120 -

2


Numerické metody

Pˇ r´ıklad 6.2.2. Pomoc´ı sloˇzených integraˇcn´ıch pravidel (6.2.2), (6.2.3) a (6.2.4) vypoˇctˇete pˇribliˇzné hodnoty integrálu

1

I=

ex dx

−1

s krokem h = 0.5, h = 0.25, h = 0.125 a h = 0.0625 a porovnejte je s pˇresnou hodnotou. ˇ sen´ı: Reˇ

Pˇresná hodnota je I = 2.350402387. Pˇribliˇzné hodnoty a jejich po-

rovnán´ı s pˇresnou hodnotou uvád´ıme v tabulce 6.2.1. Odtud je vidˇet, ˇze sloˇzené obdéln´ıkové a sloˇzené lichobˇeˇzn´ıkové pravidlo se chovaj´ı podobnˇe, zat´ımco sloˇzené Simpsonovo pravidlo dává v´ ysledky, které jsou pˇresnˇejˇs´ı o nˇekolik ˇra´d˚ u. Tabulka 6.2.1: Sloˇzené integraˇcn´ı vzorce. h

ISO

|I − ISO |

ISL

|I − ISL |

ISS

|I − ISS |

0.5 0.25 0.125 0.0625

2.326096 2.344293 2.348873 2.350020

0.024306 0.006110 0.001530 0.000383

2.399166 2.362631 2.353462 2.351167

0.048764 0.012229 0.003060 0.000765

2.351195 2.350453 2.350406 2.350402

0.000792 0.000051 0.000003 0.000000

Kontroln´ı ot´ azky Otázka 1. Jak se odvozuj´ı sloˇzené integraˇcn´ı vzorce? Otázka 2. Jakého ˇra´du jsou základn´ı sloˇzená integraˇcn´ı pravidla? ´ Ulohy k samostatn´ emu ˇreˇsen´ı 1.

Pomoc´ı sloˇzeného obdéln´ıkového, lichobˇeˇzn´ıkového a Simpsonova pravidla

vypoˇctˇete pˇribliˇzné hodnotu integrál˚ u 1 2 a) ex dx pro h = 0.1; b) 0

π

cos x2 dx pro h = π/10.

0

V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. a) ISO = 1.460393, ISL = 1.467175, ISS = 1.462681; b) ISO = 0.553752, ISL = 0.588876, ISS = 0.566030.

- 121 -


6.3.

Numerické metody

V´ ypoˇ cet integr´ alu se zadanou pˇ resnost´ı

C´ıle U integraˇcn´ıch pravidel z pˇredchoz´ıho odstavce nen´ı jasné jak dosáhnout zadané pˇresnosti. Pro tyto u ´ˇcely pouˇzijeme dvojný pˇrepoˇcet, kter´ y lze kombinovat s extrapolaˇcn´ımi vzorci.

Pˇredpokl´ adan´ e znalosti Sloˇzená integraˇcn´ı pravidla a jejich ˇra´d. Taylor˚ uv rozvoj.

V´ yklad Ve vˇetˇe 6.2.1. jsme uvedli vzorce pro odhad chyby sloˇzených integraˇcn´ıch pravidel. Jejich pouˇzit´ım lze urˇcit krok h, kter´ y zajist´ı, aby chyba vypoˇctené hodnoty integrálu byla menˇs´ı neˇz je zadaná tolerance. Tento postup má vˇsak dva háˇcky. Pˇredevˇs´ım (jak je patrné z d˚ ukazu zm´ınˇené vˇety) v´ yraz pro odhad chyby závis´ı na derivac´ıch integrované funkce, které se odhaduj´ı pracnˇe. Druh´ y problém spoˇc´ıvá v tom, ˇze výsledné odhady jsou pesimistické, takˇze integrály jsou poˇc´ıtány mnohem pˇresnˇeji neˇz je poˇzadováno, coˇz vyˇzaduje zbyteˇcnˇe velk´ y objem v´ ypoˇct˚ u. Obvykle se proto postupuje tak, ˇze hledaný integrál vyˇcislujeme opakovanˇe, stále pˇresnˇeji a ze shody v´ ysledk˚ u se usuzuje, zda jiˇz byla dosaˇzena poˇzadovaná pˇresnost. Pro tyto u ´ˇcely se osvˇedˇcilo postupné zdvojnásobován´ı poˇctu d´ılk˚ u, na nˇeˇz rozdˇelujeme integraˇcn´ı interval a, b. Hovoˇr´ıme pak o dvojném pˇrepoˇctu. Celý postup zap´ıˇseme v podobˇe algoritmu, kde I(h) oznaˇcuje pˇribliˇznou hodnotu integrálu vypoˇctenou s krokem h pomoc´ı nˇekterého sloˇzeného integraˇcn´ıho pravidla. Na zaˇcátku pˇredpokládáme, ˇze integraˇcn´ı interval a, b se rozdˇel´ı na m u ´ sek˚ u. Symbolem oznaˇcujeme zadanou pˇresnost.

- 122 -


Numerické metody

Algoritmus (Dvojn´ y pˇ repoˇ cet) Vstup: f , a, b, m, . Vypoˇcti h := (b − a)/m a I(h). Opakuj: poloˇz h := h/2, vypoˇcti I(h), dokud |I(2h) − I(h)| > . Výstup: I(h).

Pˇ r´ıklad 6.3.1. Pomoc´ı dvojného pˇrepoˇctu vypoˇctˇete pˇribliˇznou hodnotu integrálu

1

I=

ex dx

−1

s pˇresnost´ı = 0.0001. Pouˇzijte sloˇzené lichobˇeˇzn´ıkové pravidlo a zaˇcnˇete s rozdˇelen´ım integraˇcn´ıho intervalu na m = 4 d´ılk˚ u. ˇ sen´ı: V´ Reˇ ypoˇcet je zaznamenán v tabulce 6.3.1. Jako v´ ysledek dostáváme I = 2.3504 ± 0.0001.

2 Tabulka 6.3.1: Dvojn´ y pˇrepoˇcet. m

h

I(h)

|I(2h) − I(h)|

4 8 16 32 64 128 256

0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250 0.0078125

2.3991662 2.3626313 2.3534620 2.3511674 2.3505936 2.3504502 2.3504143

— 0.0365349 0.0091693 0.0022945 0.0005737 0.0001434 0.0000358

V dalˇs´ım ukáˇzeme efektivnˇejˇs´ı variantu dvojného pˇrepoˇctu s extrapolaˇcn´ımi vzorci, které nejdˇr´ıve odvod´ıme. Na hodnotu I(h) m˚ uˇzeme pohl´ıˇzet jako na funkci promˇenné h, takˇze pro ni m˚ uˇzeme psát Taylor˚ uv rozvoj. Odhady pro chybu podle vˇety 6.2.1. ukazuj´ı, ˇze se v tomto rozvoji budou vyskytovat mocniny h

- 123 -


Numerické metody

odpov´ıdaj´ıc´ı ˇra´du pˇr´ısluˇsného integraˇcn´ıho pravidla a vyˇsˇs´ı, tj. I(h) = I + Chp + O(hr ),

(6.3.1)

kde p je ˇra´d a r > p. Jestliˇze vzorec (6.3.1) nap´ıˇseme pro 2h dostaneme I(2h) = I + 2p Chp + O(hr ),

(6.3.2)

protoˇze O((2h)r ) = O(hr ). Vylouˇc´ıme-li z rovnost´ı (6.3.1) a (6.3.2) v´ yraz Chp dostáváme I = I(h) +

I(h) − I(2h) + O(hr ). 2p − 1

(6.3.3)

Odtud vid´ıme následuj´ıc´ı: 1) Veliˇcina I1 (h) = I(h) +

I(h) − I(2h) 2p − 1

(6.3.4)

je lepˇs´ı aproximac´ı I neˇz I(h) (je vyˇsˇs´ıho ˇra´du r). 2) V´ yraz E(h) =

I(h) − I(2h) 2p − 1

(6.3.5)

je aproximace chyby pˇribliˇzné hodnoty integrálu I(h), kterou m˚ uˇzeme vzhledem k bodu 1) pouˇz´ıt také pro odhad chyby aproximace I1 (h). Postup pro zvyˇsován´ı pˇresnosti v´ ysledk˚ u z´ıskan´ ych numerickou metodou podle vzorce (6.3.3) se naz´ yvá Richardsonova extrapolace. Pomoc´ı extrapolaˇcn´ıch vzorc˚ u (6.3.3) a (6.3.4) uprav´ıme algoritmus dvojného pˇrepoˇctu. Algoritmus (Dvojn´ y pˇ repoˇ cet s extrapolac´ı) Vstup: f , a, b, m, . Vypoˇcti h := (b − a)/m a I(h). Opakuj: poloˇz h := h/2, vypoˇcti I(h), dokud |E(h)| > . Výstup: I1 (h).

- 124 -


Numerické metody

Pozn´ amka V algoritmu potˇrebujeme zn´ at ˇra´d p pouˇz´ıvaného integraˇcn´ıho vzorce. Podle vˇety 6.2.1. je p = 2 pro sloˇzené obdéln´ıkové a lichobˇeˇzn´ıkové pravidlo a p = 4 pro sloˇzené Simpsonovo pravidlo. Pˇ r´ıklad 6.3.2. Pomoc´ı dvojného pˇrepoˇctu s extrapolac´ı vypoˇctˇete pˇribliˇznou hodnotu integrálu

1

I=

ex dx

−1 −4

s pˇresnost´ı = 10 . Pouˇzijete sloˇzené lichobˇeˇzn´ıkové pravidlo a zaˇcnˇete pro m = 4. ˇ sen´ı: Protoˇze p = 2, budou m´ıt extrapolaˇcn´ı vzorce tvar Reˇ I1 (h) = I(h) +

I(h) − I(2h) 3

a E(h) =

I(h) − I(2h) . 3

Pr˚ ubˇeh v´ ypoˇctu je zaznamenán v tabulce 6.3.2. Tabulka 6.3.2: Dvojn´ y pˇrepoˇcet s extrapolac´ı, lichobˇeˇzn´ıkové pravidlo. m

h

I(h)

|E(h)|

4 8 16 32 64 128

0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250

2.3991662 2.3626313 2.3534620 2.3511674 2.3505936 2.3504502

— 0.0121783 0.0030564 0.0007649 0.0001913 0.0000478

Z posledn´ıch dvou hodnot ve sloupci I(h) vypoˇc´ıtáme zpˇresnˇenou aproximaci I1 (h) = 2.3504502 +

2.3504502 − 2.3505936 = 2.3504024. 3

Jako v´ ysledek dostáváme I = 2.3504024 ± 0.0001. Porovnán´ım s pˇresnou hodnotou integrálu I = 2.350402387... vid´ıme, ˇze dosaˇzená pˇresnost je velmi vysoká. Pˇ r´ıklad 6.3.3. V pˇredchoz´ım pˇr´ıkladu pouˇzijte sloˇzené Simpsonovo pravidlo.

- 125 -


Numerické metody

ˇ sen´ı: Pouˇzijeme extrapolaˇcn´ı vzorce Reˇ I1 (h) = I(h) +

I(h) − I(2h) 15

a E(h) =

I(h) − I(2h) . 15

Pr˚ ubˇeh v´ ypoˇctu je zaznamenán v tabulce 6.3.3. Tabulka 6.3.3: Dvojn´ y pˇrepoˇcet s extrapolac´ı, Simpsonovo pravidlo. m 4 8

h

|E(h)|

I(h)

0.5000000 2.3511948 — 0.2500000 2.3504530 0.000049

Pomoc´ı extrapolace vypoˇc´ıtáme zpˇresnˇenou aproximaci I1 (h) = 2.3504530 +

2.3504530 − 2.3511948 = 2.350403. 15

Porovnán´ım s v´ ypoˇctem v pˇredchoz´ım pˇr´ıkladu vid´ıme, ˇze poˇzadované pˇresnosti jsme dosáhli mnohem dˇr´ıve. Kontroln´ı ot´ azky Otázka 1. Vysvˇetlete princip dvojného pˇrepoˇctu. Otázka 2. Jak vzniknou extrapolaˇcn´ı vzorce a jakou informaci poskytuj´ı?

´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Pomoc´ı dvojného pˇrepoˇctu vypoˇctˇete hodnoty integrál˚ u 1 π 2 a) ex dx; b) cos x2 dx 0

0

s pˇresnost´ı = 0.0001. Pouˇzijte sloˇzené lichobˇeˇzn´ıkové pravidlo. 2. V pˇredchoz´ı u ´ loze pouˇzijte sloˇzené Simpsonovo pravidlo a extrapolaci. V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. Poˇzadovanou pˇresnost dosáhneme: a) pro m = 128, kdy I = 1.462679 ± 10−4 ; b) pro m = 512, kdy I = 0.565702 ± 10−4. 2. Poˇzadovanou pˇresnost dosáhneme: a) pro m = 8, kdy I = 1.462658 ± 10−4 ; b) pro m = 32, kdy I = 0.565689 ± 10−4 .

- 126 -

Numerické metody

6.4.


Numerick´ e derivov´ an´ı

C´ıle ych Ukáˇzeme jak pˇribliˇznˇe vypoˇc´ıtat hodnoty derivac´ı f (x) a f (x) ze znám´ funkˇcn´ıch hodnot f (x − h), f (x) a f (x + h).

Pˇredpokl´ adan´ e znalosti Newton˚ uv tvar interpolaˇcn´ıho polynomu. Taylor˚ uv rozvoj.

V´ yklad Budeme postupovat podobnˇe jako pˇri odvozen´ı Newton-Cotesov´ ych vzorc˚ u v odstavci 6.1. K funkci f sestav´ıme interpolaˇcn´ı polynom pn a ten pak derivujeme m´ısto f . Pouˇzijeme pˇritom Newton˚ uv tvar interpolaˇcn´ıho polynomu (5.1.6). 1) Pro uzly x0 = x a x1 = x + h sestav´ıme lineárn´ı interpolaˇcn´ı polynom a vyjádˇr´ıme jeho prvn´ı derivaci: p1 (t) = f (x) + f [x + h, x](t − x)

⇒

p1 (t) = f [x + h, x].

Z definice pomˇern´ ych diferenc´ı dostaneme p1 (x) =

f (x + h) − f (x) ≈ f (x). h

(6.4.1)

y interpolaˇcn´ı 2) Pro uzly x0 = x−h, x1 = x a x2 = x+h sestav´ıme kvadratick´ polynom: p2 (t) = f (x − h) + f [x, x − h](t − x + h) + f [x + h, x, x − h](t − x + h)(t − x). Prvn´ı a druhá derivace maj´ı tvar p2 (t) = f [x, x − h] + f [x + h, x, x − h](2t − 2x + h), p2 (t) = 2f [x + h, x, x − h].

- 127 -


Numerické metody

Dosazen´ım t = x a u ´ pravou dostáváme p2 (x) = =

f (x) − f (x − h) f (x + h) − 2f (x) + f (x − h) + h h 2h2 f (x + h) − f (x − h) ≈ f (x) 2h

(6.4.2)

a p2 (x) =

f (x + h) − 2f (x) + f (x − h) ≈ f (x). 2 h

(6.4.3)

Pˇ r´ıklad 6.4.1. Vypoˇctˇete pˇribliˇzné hodnoty prvn´ı a druhé derivace pro funkci f (x) = sin x v bodˇe x = 1 s krokem h = 0.01 pomoc´ı vzorc˚ u (6.4.1), (6.4.2) a (6.4.3). Porovnejte je s pˇresnými hodnotami. ˇ sen´ı: Reˇ

Budeme potˇrebovat funkˇcn´ı hodnoty sin 1 = 0.841471, sin 1.01 =

0.846832 a sin 0.99 = 0.836026. Podle vzorec (6.4.1), resp. (6.4.2) dostaneme f (1) ≈ p1 (1) =

0.846832 − 0.841471 = 0.536100, 0.01

f (1) ≈ p2 (1) =

0.846832 − 0.836026 = 0.540300. 2 · 0.01

resp.

Pˇresná hodnota je f (1) = cos 1 = 0.540302, takˇze druhá pˇribliˇzná hodnota je podstatnˇe pˇresnˇejˇs´ı. Podle vzorce (6.4.3) vypoˇc´ıtáme druhou derivaci f (1) ≈ p2 (1) =

0.846832 − 2 · 0.841471 + 0.836026 = −0.840000. 0.012

Nyn´ı je pˇresná hodnota f (1) = − sin 1 = −0.841471. Následuj´ıc´ı vˇeta ukazuje jak´ y je ˇra´d jednotliv´ ych vzorc˚ u.

- 128 -

Numerické metody


Vˇ eta 6.4.1. Necht’ je funkce f dostateˇcnˇe hladká (má v okol´ı bodu x spojitou druhou, tˇret´ı, resp. ˇctvrtou derivaci). Pak plat´ı:

f (ξ) , 2 f (3) (ξ) , p2 (x) − f (x) = h2 3 f (4) (ξ) p2 (x) − f (x) = h2 . 12 p1 (x) − f (x) = h

(6.4.4)

Ve vˇsech pˇr´ıpadech je ξ bl´ıˇze neurˇcený bod z okol´ı x. D˚ ukaz: Nejjednoduˇsˇs´ı zp˚ usob jak z´ıskat tato tvrzen´ı je pouˇz´ıt Taylor˚ uv rozvoj. Napˇr´ıklad pro odvozen´ı (6.4.4) staˇc´ı vyjádˇrit prvn´ı derivaci z

2f

f (x + h) = f (x) + hf (x) + h

(ξ) . 2 2

Výpoˇcet pˇribliˇzných hodnot derivac´ı mohou podstatnˇe ovlivnit zaokrouhlovac´ı chyby. Jmenovatele vzorc˚ u totiˇz obsahuj´ı parametr h, kter´ y mus´ı b´ yt mal´ y, abychom dostali dostateˇcnˇe pˇresnou aproximaci derivace. Souˇcasnˇe ovˇsem malá hodnota jmenovatele zlomku zvˇetˇsuje zaokrouhlovac´ı chyby v ˇcitateli. Nejvˇetˇs´ı dosaˇzitelná pˇresnost proto mus´ı b´ yt jist´ ym kompromisem mezi chybou aproximace a chybou zaokrouhlen´ı. Na pˇr´ıkladu vzorce (6.4.1) si ukáˇzeme anal´ yzu tohoto problému. y vznikne jako souˇcet odOznaˇcme Ecelk (h) horn´ı odhad celkové chyby, kter´ hadu chyby aproximace e(h) a chyby zaokrouhlen´ı z(h), tj. Ecelk (h) = e(h) + z(h). Podle (6.4.4) je e(h) = Ch. Dále budeme pˇredpokládat, ˇze pˇri v´ ypoˇctu f (x) dostaneme vlivem zaokrouhlen´ı poruˇsenou hodnotu f ∗ (x), pˇriˇcemˇz velikost poruchy je nejvýˇse κ, tj. |f ∗ (x) − f (x)| ≤ κ. Pro vzorec (6.4.1) dostáváme následuj´ıc´ı odhad:

- 129 -


Numerické metody

f (x + h) − f (x) f ∗ (x + h) − f ∗ (x)

≤

−

h h ≤

1 2κ (|f ∗(x) − f (x)| + |f ∗ (x + h) − f (x + h)|) ≤ = z(h). h h

Odhad celkové chyby má proto tvar Ecelk (h) = Ch+ 2κ , viz obrázek 6.4.1. Funkce h Ecelk (h) má jediné minimum v bodˇe ' hopt =

2κ C

(6.4.5)

a odpov´ıdaj´ıc´ı hodnota Ecelk (hopt ) pˇredstavuje nejmenˇs´ı horn´ı odhad celkové chyby. Pˇri v´ ypoˇctu derivace podle (6.4.1) tedy dojde pro h menˇs´ı neˇz hopt paradoxnˇe ke ztrátˇe pˇresnosti.

Ecelk (h) e(h)

z(h) hopt

h

Obrázek 6.4.1: Odhad celkové chyby Ecelk (h).

Pˇ r´ıklad 6.4.2. Vypoˇctˇete pˇribliˇznou hodnotu prvn´ı derivace funkce f (x) = sin x v bodˇe x = 1 s krokem h = 0.01 a h = 0.001 pomoc´ı vzorce (6.4.1). Zaokrouhlujte pˇritom na ˇctyˇri desetinná m´ısta. V´ ysledky porovnejte s optimáln´ım krokem hopt .

- 130 -


Numerické metody ˇ sen´ı: Reˇ

V naˇsem pˇr´ıpadˇe je κ = 5 · 10−5 a

f (ξ) − sin ξ

=

2 2 ≤ 0.5 = C.

Dosazen´ım do vzorce (6.4.5) vypoˇc´ıtáme hopt = 0.014142, takˇze výsledek pro h = 0.01 by mˇel být pˇresnˇejˇs´ı. Pˇribliˇzné derivace maj´ı hodnotu: f (1) ≈

sin(1.01) − sin(1) 0.01

f (1) ≈

sin(1.001) − sin(1) . 0.8420 − 0.8415 = = 0.50. 0.001 0.001

. 0.8468 − 0.8415 = 0.53, = 0.01

Jejich porovnán´ı s hodnotou f (1) = cos 1 = 0.540302 potvrzuje pˇredchoz´ı závˇer.

Kontroln´ı ot´ azky Otázka 1. Jak se odvozuj´ı vzorce pro pˇribliˇzný v´ ypoˇcet derivace? Otázka 2. Znázornˇete graficky smysl vzorc˚ u pro v´ ypoˇcet prvn´ı derivace. Otázka 3. Jak ovlivˇ nuj´ı v´ ypoˇcet derivac´ı zaokrouhlovac´ı chyby? Otázka 4. Odvod’te podrobnˇe vztah (6.4.5). ´ Ulohy k samostatn´ emu ˇreˇsen´ı 1.

Vypoˇctˇete pˇribliˇznˇe prvn´ı a druhou derivaci funkce f (x) = cos x v bodˇe

x = 1.5 s krokem h = 0.01. V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. f (1.5) ≈ −0.997832 podle vzorce (6.4.1); f (1.5) ≈ −0.997478 podle vzorce (6.4.2); f (1.5) ≈ −0.070737 podle vzorce (6.4.3). Shrnut´ı lekce Ukázali jsme pouˇzit´ı základn´ıch integraˇcn´ıch pravidel a vzorc˚ u numerického derivován´ı. Vidˇeli jsme, ˇze numerické v´ ypoˇcty integrál˚ u jsou stabiln´ı vzhledem k zaokrouhlovac´ım chybám. Zat´ımco pˇribliˇzné v´ ypoˇcty derivac´ı jsou na zaokrouhlovac´ı chyby pomˇernˇe citlivé.

- 131 -

integrály lze vypočítat snadno pomocí tabulek a klasických integračních metod jako je per partes nebo substituce. Tak například integrály

Recommend Documents