ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
7.
Numerick´e metody
ODR – poˇ c´ ateˇ cn´ı u ´ lohy
Pr˚ uvodce studiem Jen velmi m´alo diferenci´aln´ıch rovnic, kter´e se vyskytuj´ı pˇri popisu praktick´ ych u ´ loh, se d´a ˇreˇsit exaktnˇe, a i kdyˇz dok´aˇzeme naj´ıt vzorce popisuj´ıc´ı analytick´e ˇreˇsen´ı, b´ yvaj´ı natolik sloˇzit´e, ˇze je obt´ıˇzn´e s nimi d´ale pracovat. Numerick´e metody jsou proto ˇcasto jedinou moˇznost´ı jak danou diferenci´aln´ı rovnici vyˇreˇsit. V t´eto kapitole pˇredstav´ıme z´akladn´ı numerick´e metody pro ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh pro obyˇcejn´e diferenci´aln´ı rovnice (ODR). Nejdˇr´ıve si na pˇr´ıkladu pˇripomeneme nˇekter´e analytick´e postupy a nachyst´ame si modelovou u ´ lohu. Pak struˇcnˇe odvod´ıme nejjednoduˇsˇs´ı Eulerovu metodu, na n´ıˇz uk´aˇzeme podstatn´e rysy spoleˇcn´e vˇsem numerick´ ym metod´am. U ostatn´ıch metod pak uvedeme pouze vzorce a na pˇr´ıkladech pˇredvedeme jejich pouˇzit´ı a vlastnosti.
7.1.
Formulace u ´lohy
C´ıle Zavedeme obecnou formulaci poˇc´ateˇcn´ı u ´ lohy pro ODR prvn´ıho ˇra´du a zm´ın´ıme pˇredpoklady, kter´e zajiˇst’uj´ı existenci a jednoznaˇcnost ˇreˇsen´ı. Na pˇr´ıkladu pˇripomeneme exaktn´ı metody pro v´ ypoˇcet analytick´eho ˇreˇsen´ı.
Pˇredpokl´ adan´ e znalosti Exaktn´ı metody ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh pro ODR.
- 132 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
V´ yklad Budeme se zab´yvat numerick´ ym v´ ypoˇctem funkce y = y(x), kter´a na intervalu a, b vyhovuje rovnici y (x) = f (x, y(x)),
(7.1.1)
kde f = f (x, y) je dan´a prav´a strana. Rovnice (7.1.1) je obyˇcejn´ a diferenci´ aln´ı rovnice prvn´ıho ˇra´du a z teorie je o n´ı zn´amo, ˇze ˇreˇsen´ı (pokud existuje) je urˇceno aˇz na jeden voliteln´y parametr. Aby bylo ˇreˇsen´ı urˇceno jednoznaˇcnˇe, budeme poˇzadovat splnˇen´ı poˇc´ateˇcn´ı podm´ınky ve tvaru y(a) = c.
(7.1.2)
´ Uloha (7.1.1), (7.1.2) je poˇca´teˇcn´ı u ´loha pro obyˇcejnou diferenci´aln´ı rovnici. Ovˇsem ani tato u ´ loha nemus´ı m´ıt jedin´e ˇreˇsen´ı. Vˇse z´avis´ı na vlastnostech prav´e strany f . O t´eto funkci budeme v cel´e kapitole pˇredpokl´adat, ˇze je spojit´a v prvn´ı nuje Lipschitzovu podm´ınku ve druh´e promˇenn´e, tj. existuje konpromˇenn´e a splˇ stanta L (nez´avisl´a na x a y) takov´a, ˇze |f (x, y) − f (x, z)| ≤ L|y − z| ∀x ∈ a, b ∀y, z ∈ R.
(7.1.3)
Za tˇechto pˇredpoklad˚ u m´a poˇc´ateˇcn´ı u ´ loha jedin´e ˇreˇsen´ı. Pˇ r´ıklad 7.1.1. Poˇc´ateˇcn´ı u ´ loha y = y2,
y(0) = 1
nem´a jedin´e ˇreˇsen´ı na ntervalu 0, 2, protoˇze prav´a strana na uveden´em intervalu nesplˇ nuje Lipschitzovu podm´ınku. V n´asleduj´ıc´ım pˇr´ıkladu si pˇripomeneme klasick´e postupy pro ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh. Bude se jednat o u ´ lohu, kterou v dalˇs´ım textu pouˇzijeme jako modelov´y pˇr´ıklad pˇri testov´an´ı numerick´ ych metod.
- 133 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Pˇ r´ıklad 7.1.2.
Numerick´e metody
Ovˇeˇrte, ˇze poˇc´ateˇcn´ı u ´ loha y = x2 − 0.2y,
y(−2) = −1
(7.1.4)
m´a jedin´e ˇreˇsen´ı. Toto ˇreˇsen´ı vypoˇctˇete exaktnˇe pomoc´ı zn´am´ ych analytick´ ych metod. ˇ sen´ı: Reˇ
Ovˇeˇren´ı Lipschitzovy podm´ınky pro f (x, y) = x2 − 0.2y provedeme
dosazen´ım do x2 − 0.2y − x2 + 0.2z f (x, y) − f (x, z) = = −0.2, y−z y−z ˇ sen´ı odkud vid´ıme, ˇze (7.1.3) je splnˇeno napˇr´ıklad pro konstantu L = 0.2. Reˇ zadan´e poˇc´ateˇcn´ı u ´ lohy je proto jedin´e. Pˇri jeho v´ ypoˇctu nejdˇr´ıve vyˇreˇs´ıme homogenn´ı diferenci´aln´ı rovnici u = −0.2u pomoc´ı separace promˇenn´ych. Postup je n´asleduj´ıc´ı: du du = −0.2u ⇒ = −0.2 dx ⇒ dx u
du = u
−0.2 dx ⇒
ln |u| = −0.2x + c ⇒ |u| = e−0.2x+c ⇒ u(x) = Ce−0.2x ,
C > 0.
Pro v´ ypoˇcet (partikul´arn´ıho) ˇreˇsen´ı yp rovnice (7.1.4) pouˇzijeme metodu variace konstant, kdy pˇredem pˇredpokl´ad´ame, ˇze ˇreˇsen´ı bude ve tvaru yp (x) = C(x)e−0.2x . Dosazen´ım do diferenci´aln´ı rovnice (7.1.4) dostaneme C (x)e−0.2x + C(x)(−0.2)e−0.2x = x2 − 0.2C(x)e−0.2x ⇒ C (x) = x2 e0.2x a odtud integrac´ı per partes vypoˇcteme C(x) = x2 e0.2x dx = e0.2x (5x2 − 50x + 250), takˇze yp (x) = 5x2 −50x+ 250. Obecn´e ˇreˇsen´ı dan´e diferenci´aln´ı rovnice lze zapsat ve tvaru y(x) = yp (x) + u(x), tj. y(x) = 5x2 − 50x + 250 + Ce−0.2x . Koneˇcnˇe
- 134 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody z poˇc´ateˇcn´ı podm´ınky urˇc´ıme konstantu C
−1 = y(−2) = 20 + 100 + 250 + Ce0.4 ⇒ C = −248.688737 a t´ım i ˇreˇsen´ı naˇs´ı v´ ychoz´ı u ´ lohy: y(x) = 5x2 − 50x + 250 − 248.018417e−0.2x .
(7.1.5)
Pozn´ amka ´ Ulohu (7.1.1), (7.1.2) m˚ uˇzeme ch´ apat vektorovˇe: hled´ ame vektorovou funkci a strana je rey(x) = (y1 (x), . . . , yn (x)) vyhovuj´ıc´ı rovnici (7.1.1), kde prav´ nuj´ıc´ı prezentov´ ana vektorovu funkc´ı f (x, y) = (f1 (x, y), . . . , fn (x, y)) , a splˇ poˇca´teˇcn´ı podm´ınku (7.1.2) zadanou vektorem c = (c1 , . . . , cn ) . Vˇsechny mumerick´e metody, s nimiˇz se v dalˇs´ım sezn´am´ıme, lze proto pouˇz´ıt i pro soustavy obyˇcejn´ych diferenci´ aln´ıch rovnic prvn´ıho ˇra´du.
Pozn´ amka Skuteˇcnost, ˇze se omezujeme pouze na rovnice prvn´ıho ˇra´du, nen´ı ˇz´adn´e podstatn´e omezen´ı. Z teorie diferenci´ aln´ıch rovnic je totiˇz zn´ amo, ˇze poˇc´ateˇcn´ı u ´lohu pro rovnici vyˇsˇs´ıho ˇra´du neˇz 1 lze jednoduchou substituc´ı pˇrev´est na soustavu rovnic prvn´ıho ˇra´du.
Kontroln´ı ot´ azky Ot´azka 1. Jak vypad´a obecn´a formulace poˇc´ateˇcn´ı u ´ lohy pro obyˇcejnou diferenci´aln´ı rovnici prvn´ıho ˇra´du? Ot´azka 2. Zopakujte si exaktn´ı metody pro v´ ypoˇcet pˇresn´eho (analytick´eho) ˇreˇsen´ı.
- 135 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
7.2.
Numerick´e metody
Eulerova metoda
C´ıle Odvod´ıme rekurentn´ı vzorce pro Eulerovu metodu a uk´aˇzeme jejich pouˇzit´ı.
Pˇredpokl´ adan´ e znalosti Prov´adˇen´ı rekurentn´ıch v´ ypoˇct˚ u, analytick´e ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh pro ODR.
V´ yklad Z´akladem, z nˇehoˇz vych´azej´ı numerick´e metody je diskretizace promˇenn´ych. Pˇribliˇzn´e ˇreˇsen´ı se nekonstruuje jako spojit´a funkce, ale postupnˇe se generuje posloupnost uzl˚ u x0 = a, x1 , x2 , . . . a pro nˇe se stanov´ı ˇc´ısla y0 = c, y1 , y2 , . . ., kter´a aproximuj´ı hodnoty pˇresn´eho ˇreˇsen´ı y(x0 ), y(x1 ), y(x2 ), . . .. Pro jednoduchost budeme pˇredpokl´adat, ˇze uzly jsou ekvidistantn´ı s krokem h a plat´ı xn = b, tj. h = (b − a)/n a xi = a + ih, i = 0, 1, . . . , n. Pˇri odvozov´an´ı vzorc˚ u pro pˇribliˇzn´e ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh m˚ uˇzeme postupovat napˇr´ıklad tak, ˇze v diferenci´aln´ı rovnici (7.1.1) zapsan´e v uzlu xi nahrad´ıme pˇresnou hodnotu y(xi) jej´ı aproximac´ı yi a derivaci na lev´e stranˇe vyj´adˇr´ıme numerick´ ym vzorcem (6.4.1), tj. y (xi ) = f (xi , y(xi))
≈
yi+1 − yi = f (xi , yi). h
Odtud jsme schopni pˇri zn´am´ ych hodnot´ach xi a yi vypoˇc´ıtat yi+1 . Po pˇrepisu do explicitn´ıho tvaru dost´av´ame rekurentn´ı vzorce Eulerovy metody: ⎫ ⎬ y0 = c, = y + hf (x , y ), i = 0, 1, . . . , n − 1. ⎭ y i+1
i
i
(7.2.1)
i
Pˇ r´ıklad 7.2.1. Poˇc´ateˇcn´ı u ´ lohu (7.1.4) ˇreˇste na intervalu −2, 3 pomoc´ı Eulerovy metody s krokem h = 1 a 0.5. V´ ysledky porovnejte s pˇresn´ym ˇreˇsen´ım.
- 136 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
ˇ sen´ı: V naˇsem pˇr´ıpadˇe je x0 = −2, y(−2) = y0 = −1 a f (x, y) = x2 − 0.2y. Reˇ Pro h = 1 dost´av´ame: y1 = −1 + 1 · ((−2)2 − 0.2 · (−1)) = 3.2, y2 = 3.2 + 1 · ((−1)2 − 0.2 · 3.2) = 3.56, y3 = 2.848, y4 = 3.2784, y5 = 6.6227. Podobnˇe pro h = 0.5 vypoˇcteme: y1 = −1 + 0.5 · ((−2)2 − 0.2 · (−1)) = 1.1, y2 = 1.1 + 0.5 · ((−1.5)2 − 0.2 · 1.1) = 2.115, y3 = 2.4035, . . . , y10 = 7.4988. Pro lepˇs´ı n´azornost zapiˇsmˇe obˇe pˇribliˇzn´a ˇreˇsen´ı do tabulky 7.2.1, v n´ıˇz provedeme tak´e porovn´an´ı s pˇresn´ym ˇreˇsen´ım y(x) urˇcen´ym vzorcem (7.1.5). Z tabulky je Tabulka 7.2.1: Eulerova metoda. h = 0.5
h=1 i
xi
yi
|yi − y(xi)|
0
−2
−1
0
1
−1
3.2
1.9491
2
0
3.56
2.2487
3
1
2.848
1.4571
4
2
3.2784
0.0206
5
3
6.6227
1.8940
i
xi
yi
0 1 2 3 4 5 6 7 8 9 10
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3
−1 1.1 2.1150 2.4035 2.2882 2.0593 1.9784 2.2806 3.1775 4.8598 7.4988
|yi − y(xi )| 0 0.5447 0.8641 0.9971 0.9769 0.8322 0.5875 0.2637 0.1214 0.5529 1.0179
vidˇet, ˇze ˇreˇsen´ı vypoˇcten´e s menˇs´ım krokem h je pˇresnˇejˇs´ı. Tento z´avˇer potvrzuje tak´e obr´azek 7.2.1, kde jsou pˇribliˇzn´a ˇreˇsen´ı zn´azornˇena jako lomen´e ˇc´ary.
- 137 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
9 8 7 6 5 4
h=1
3
h = 0.5
2 1
y(x) 0 −1 −2
−1
0
1
2
3
Obr´azek 7.2.1: Pˇribliˇzn´a ˇreˇsen´ı z Eulerovy metody a pˇresn´e ˇreˇsen´ı (teˇckovanˇe). Pozn´ amka Eulerova metoda je metodou prvn´ıho ˇra´du, pˇrestoˇze jsme ji odvodili ze vzorec (6.4.1), kter´y je ˇra´du druh´eho. Pˇri v´ypoˇctu podle Eulerovy metody doch´ az´ı ke kumuluaci (diskretizaˇcn´ı) chyby, coˇz m´ a za d˚ usledek sn´ıˇzen´ı ˇra´du. Kontroln´ı ot´ azky Ot´azka 1. Odvod’te Eulerovu metodu. Jak´eho je ˇra´du a proˇc? ´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Eulerovou metodou ˇreˇste diferenci´aln´ı rovnici y = 1/(x2 +1)−0.1y s poˇc´ateˇcn´ı podm´ınkou y(−2) = −1 na intervalu −2, 3 s krokem h = 0.5. V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. y0 = −1, y1 = −0.85, y2 = −0.6537, y3 = −0.3710, y4 = 0.0476, y5 = 0.5452, y6 = 0.9179, y7 = 1.1220, y8 = 1.2198, y9 = 1.2588, y10 = 1.2648.
- 138 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
7.3.
Jednokrokov´ e metody vyˇ sˇ s´ıho ˇ r´ adu
C´ıle Uk´aˇzeme pouˇzit´ı rekurentn´ıch vzorc˚ u Heunovy metody a klasick´e RungeovyKuttovy metody.
Pˇredpokl´ adan´ e znalosti Eulerova metoda, analytick´e ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh pro ODR.
V´ yklad N´asleduj´ıc´ı vzorce popisuj´ı metody vyˇsˇs´ıho ˇra´du. Pˇri jejich pouˇzit´ı dos´ahneme srovnateln´e pˇresnosti jako u Eulerovy metody pˇri podstatnˇe vˇetˇs´ım kroku h, tedy pˇri v´ yraznˇe menˇs´ım celkov´em objemu v´ ypoˇct˚ u. Heunova metoda, kter´a je druh´eho ˇra´du, je urˇcena rekurentn´ımi vzorci: ⎫ ⎪ ⎪ y0 = c, ⎪ ⎪ ⎪ ⎪ ⎬ k1 = hf (xi , yi ), (7.3.1) ⎪ ⎪ k2 = hf (xi + h, yi + k1 ), ⎪ ⎪ ⎪ ⎪ 1 = y + (k + k ), i = 0, . . . , n − 1. ⎭ y i+1
i
2
1
2
Rungeova-Kuttova metoda ˇctvrt´eho ˇra´du (RK4) je urˇcena rekurentn´ımi vzorci: ⎫ ⎪ ⎪ y0 = c, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ k1 = hf (xi , yi), ⎪ ⎪ ⎪ ⎪ 1 1 ⎬ k2 = hf (xi + 2 h, yi + 2 k1 ), (7.3.2) ⎪ ⎪ k3 = hf (xi + 12 h, yi + 12 k2 ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ k4 = hf (xi + h, yi + k3 ), ⎪ ⎪ ⎪ ⎪ 1 yi+1 = yi + 6 (k1 + 2k2 + 2k3 + k4 ), i = 0, . . . , n − 1. ⎭
- 139 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
Pozn´ amka Obˇe uveden´e metody jsou jednokrokov´e, ˇc´ımˇz m´ ame namysli, ˇze pro v´ypoˇcet yi+1 staˇc´ı hodnota z jedin´eho pˇredchoz´ıho kroku yi . Pozn´ amka Ze vzorc˚ u (7.3.1), (7.3.2) je d´ ale vidˇet,ˇze zv´yˇsen´ı ˇra´du u jednokrokov´e metody je spojeno s vˇetˇs´ım objemem v´ypoˇct˚ u, zejm´ena je v kaˇzd´em kroku potˇreba poˇc´ıtat hodnoty prav´e strany f (x, y) ve v´ıce bodech.
Pˇ r´ıklad 7.3.1. Poˇc´ateˇcn´ı u ´ lohu (7.1.4) ˇreˇste na intervalu −2, 3 pomoc´ı Heunovy metody a metody RK4 s krokem h = 1. V´ ysledky porovnejte s pˇresn´ym ˇreˇsen´ım a s ˇreˇsen´ım z´ıskan´ ym Eulerovou metodou. ˇ sen´ı: Opˇet vyjdeme z x0 = −2, y(−2) = y0 = −1 a f (x, y) = x2 − 0.2y. Reˇ Naznaˇc´ıme prvn´ı dva kroky Heunovy metody: k1 = 1 · ((−2)2 − 0.2 · (−1)) = 4.2, k2 = 1 · ((−2 + 1)2 − 0.2 · (−1 + 4.2)) = 0.36, y1 = −1 + 12 (4.2 + 0.36) = 1.28, k1 = 1 · ((−1)2 − 0.2 · 1.28) = 0.744, k2 = 1 · ((−1 + 1)2 − 0.2 · (1.28 + 0.744)) = −0.4048, y2 = 1.28 + 12 (0.744 + (−0.4048)) = 1.4496, atd. Obˇe pˇribliˇzn´a ˇreˇsen´ı zap´ıˇseme do tabulky 7.3.1 a porovn´ame je s pˇresn´ym ˇreˇsen´ım y(x) urˇcen´ym vzorcem (7.1.5). Porovn´an´ım se sloupcem h = 1 v tabulce 7.2.1 vid´ıme, ˇze dosaˇzen´e v´ ysledky jsou pˇresnˇejˇs´ı neˇz u Eulerovy metody, pˇriˇcemˇz nejpˇresnˇejˇs´ı v´ ysledky d´av´a metoda RK4, viz tak´e obr´azek 7.3.1.
- 140 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
Tabulka 7.3.1: Metody vyˇsˇs´ıho ˇra´du. Metoda RK4
Heunova metoda i
xi
yi
0 −2 1 −1 2 0 3 1 4 2 5 3
−1 1.2800 1.4496 1.6887 3.7847 9.2035
|yi − y(xi )| 0 0.0291 0.1383 0.2978 0.4858 0.6867
i
xi
yi
0 −2 1 −1 2 0 3 1 4 2 5 3
−1 1.2508 1.3112 1.3910 3.2994 8.5175
|yi − y(xi )| 0 0.0001 0.0001 0.0001 0.0004 0.0008
10
8
6
4
2
y(x)
0
−2 −2
−1
0
1
2
3
Obr´azek 7.3.1: Pˇribliˇzn´e ˇreˇsen´ı vypoˇc´ıtan´e metodou RK4 (plnˇe), Heunovou metodou (ˇc´arkovanˇe) a pˇresn´e ˇreˇsen´ı (teˇckovanˇe).
- 141 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
Kontroln´ı ot´ azky Ot´azka 1. Jak´e zn´ate metody vyˇsˇs´ıho ˇra´du? Co je jejich v´ yhodou oproti metod´am prvn´ıho ˇra´du? Ot´azka 2. Jak´ y je vztah mezi ˇra´dem metody a poˇcetn´ı n´aroˇcnost´ı jednoho kroku?
´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Diferenci´aln´ı rovnici y = y(1 + sin x) − x3 , y(1) = 0 ˇreˇste na intervalu 1, 2 Heunovou metodou a metodou RK4 s krokem h = 0.2.
V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. Heunova metoda: y0 = 0, y1 = −0.3114, y2 = −0.9732, y3 = −2.2320, y4 = −4.4495, y5 = −8.1186; metoda RK4: y0 = 0, y1 = −0.3210, y2 = −1.0087, y3 = −2.3257, y4 = −4.6586, y5 = −8.5351.
- 142 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
7.4.
V´ıcekrokov´ e metody
C´ıle Uk´aˇzeme pouˇzit´ı v´ıcekrokov´ ych metod pro ˇreˇsn´ı poˇc´ateˇcn´ıch u ´ loh pro ODR. Pomoc´ı z´akladn´ıch v´ıcekrokov´ ych vzorc˚ u odvod´ıme algoritmus prediktor-korektor.
Pˇredpokl´ adan´ e znalosti Jednokrokov´e a analytick´e metody ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh pro ODR.
V´ yklad Definice 7.4.1. Dan´a metoda se na´ yv´a k-krokov´a, jestliˇze pˇri v´ ypoˇctu yi+1 je potˇreba dosadit do vzorce metody hodnoty ˇreˇsen´ı z k pˇredchoz´ıh krok˚ u, tj. hodnoty yi , yi−1 , . . . , yi−k+1. Vˇsechny dosud uveden´e metody byly jednokrokov´e. Pˇr´ıkladem dvoukrokov´e metody je vzorec yi+1 = yi−1 + 2hf (xi , yi),
(7.4.1)
kter´ y lze odvodit podobnˇe jako Eulerovu metodu (7.2.1) pomoc´ı numerick´eho derivov´an´ı. Vzorci (7.4.1) se ˇr´ık´a metoda sk´ akaj´ıc´ı ˇz´aby, protoˇze hodnoty pˇribliˇzn´eho ˇreˇsen´ı vˇetˇsinou oscilij´ı kolem ˇreˇsen´ı pˇresn´eho. Definice 7.4.2. Dan´a metoda se na´ yv´a l-bodov´a, jestliˇze pˇri v´ ypoˇctu yi+1 vyˇzaduje vzorec metody vypoˇc´ıtat hodnoty prav´e strany f v l bodech. Eulerova metoda (7.2.1) je jednobodov´a, Heunova metoda (7.3.1) je dvoubodov´a a Runge-Kuttova metoda (7.3.2) je ˇctyˇrbodov´a. Vid´ıme, ˇze u jednokro-
- 143 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
kov´ ych metod je pro dosaˇzen´ı vyˇsˇs´ıho ˇra´du poˇc´ıtat hodnoty prav´e strany ve v´ıce bodech, ˇc´ımˇz se ovˇsem zvˇetˇsuje pracnost v´ypoˇctu. Tuto nepˇr´ıjemnou vlastnost odstraˇ nuj´ı v´ıcekrokov´e metody. Hodnoty pˇribliˇzn´eho ˇreˇsen´ı z pˇredchoz´ıch krok˚ u totiˇz pˇredstavuj´ı v dan´ y okamˇzik jiˇz napoˇctenou informaci, kterou lze vyuˇz´ıt pro dosaˇzen´ı vyˇsˇs´ıho ˇra´du. Vˇsimnˇeme si, ˇze metoda (7.4.1) je jednobodov´a, dvoukrokov´a a druh´eho ˇra´du (jej´ı odvozen´ı je zaloˇzeno na pouˇzit´ı vzorce numerick´eho derivov´an´ı (6.4.2) druh´eho ˇra´du). Jistou komplikac´ı u k-krokov´e metody je skuteˇcnost, ˇze pˇri zah´ajen´ı v´ ypoˇctu ´sek. Jeho v´ ypoˇcet se potˇrebujeme zn´at hodnoty y0 , y1 , . . . , yk−1 , tzv. poˇc´ateˇcn´ı u zpravidla prov´ad´ı vhodnou jednokrokovou metodou (stejn´eho nebo vyˇsˇs´ıho ˇra´du). Vˇsechny dosud uveden´e metody byly explicitn´ı, tzn. ˇze v jejich vzorci se yi+1 vyskytovalo pouze na lev´e stranˇe. U v´ıcekrokov´ ych metod se ˇcasto pouˇz´ıvaj´ı vzorce implicitn´ı, u nichˇz se yi+1 vyskytuje i na stranˇe prav´e. Pouˇzit´ı takov´ehoto vzorce pˇredstavuje v kaˇzd´em kroku rovnici (zpravidla neline´arn´ı), jej´ımˇz ˇreˇsen´ım u vhodn´e iteraˇcn´ı metody, je yi+1. Pˇri ˇreˇsen´ı t´eto rovnice lze pouˇz´ıt nˇekolik krok˚ coˇz je podstata algoritm˚ u typu prediktor-korektor. Adamsovy-Bashforthovy vzorce Pod n´azvem Adamsovy-Bashforthovy vzorce rozumm´ıme explicitn´ı v´ıcekrokov´e metody, kter´e lze odvodit integrac´ı interpolaˇcn´ıho polynomu sestaven´eho pro derivaci ˇreˇsen´ı diferenci´aln´ı rovnice. Ukaˇzme odvozen´ı tˇr´ıkrokov´eho vzorce. Necht’ y(x) je ˇreˇsen´ı u ´ lohy (7.1.1) a oznaˇcme yi = y (xi ). Derivaci y (x) m˚ uˇzeme vyj´adˇrit ve tvaru y (x) = p2 (x) + e(x),
(7.4.2)
kde p2 (x) je Newton˚ uv tvar interpolaˇcn´ıho polynomu v uzlech xi , xi−1 , xi−2 a e(x) je interpolaˇcn´ı chyba: p2 (x) = yi + y [xi−1 , xi ](x − xi ) + y [xi−2 , xi−1 , xi ](x − xi )(x − xi−1 ),
- 144 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
e(x) =
˜ y (4) (ξ) (x − xi )(x − xi−1 )(x − xi−2 ). 3!
V (7.4.2) provedeme integraci na intervalu xi , xi+1 : xi+1 y (x) dx = y(xi+1 ) − y(xi),
xi xi+1 xi
h2 5h3 p2 (x) dx = · h + y [xi−1 , xi ] · + y [xi−2 , xi−1 , xi ] · 2 6 2 y − 2yi−1 + yi−2 5h3 y − yi−1 h · + i = yi · h + i · h 2 2h2 6 h = (23yi − 16yi−1 + 5yi−2 ), 12 yi
xi+1
e(x) dx = xi
y (4) (ξ) 9 4 3 4 (4) · · h = · h · y (ξ), 3! 4 8
takˇze dohromady dost´av´ame y(xi+1 ) − y(xi ) =
h 3 (23yi − 16yi−1 + 5yi−2 ) + h4 y (4) (ξ). 12 8
Vynech´an´ım posledn´ıho ˇclenu vznikne rekurentn´ı vzorec, kter´y lze pouˇz´ıt pro pˇribliˇzn´y v´ ypoˇcet y(xi+1 ). Staˇc´ı si uvˇedomit, ˇze (pˇribliˇzn´e) hodnoty derivac´ı , yi−2 lze snadno vypoˇc´ıtat dosazen´ım do prav´e strany f naˇs´ı diferenci´aln´ı yi , yi−1
rovnice (7.1.1). Oznaˇc´ıme-li yi ≈ y(xi) a fi = f (xi , yi) ≈ yi , obdrˇz´ıme n´asleduj´ıc´ı tvar vzorce: yi+1 = yi +
h (23fi − 16fi−1 + 5fi−2 ). 12
Jedn´a se zˇrejmˇe o tˇr´ıkrokovou metodu tˇret´ıho ˇra´du (srovnej s pozn´amkou u Eulerovy metody). Souˇcasnˇe je to metoda jednobodov´a ! Staˇc´ı si totiˇz uvˇedomit, ˇze pˇri v´ ypoˇctu yi+1 nen´ı potˇreba vyˇc´ıslit hodnoty fi−1 a fi−2 , protoˇze k tomu jiˇz doˇslo v pˇredchoz´ıch kroc´ıch. Analogick´ ym postupem lze odvodit dalˇs´ı Adamsovy-Bashforthovy vzorce liˇs´ıc´ı se poˇctem krok˚ u. Obvykle pouˇz´ıvan´e vzorce jsou nejv´yˇse ˇctyˇrkrokov´e. Zde je jejich pˇrehled: yi+1 = yi + hfi ,
(7.4.3)
- 145 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
h yi+1 = yi + (3fi − fi−1 ), 2 h yi+1 = yi + (23fi − 16fi−1 + 5fi−2 ), 12 h yi+1 = yi + (55fi − 59fi−1 + 37fi−2 − 9fi−3 ). 24
Pˇ r´ıklad 7.4.1.
(7.4.4) (7.4.5) (7.4.6)
Poˇc´ateˇcn´ı u ´ lohu (7.1.4) ˇreˇste na intervalu −2, 3 pomoc´ı
tˇr´ıkrokov´e Adamsovy-Bashforthovy metody s krokem h = 1. V´ ysledky porovnejte s pˇresn´ym ˇreˇsen´ım. ˇ sen´ı: Poˇc´ateˇcn´ı u Reˇ ´sek vypoˇc´ıt´ame metodou RK4. Jsou to prvn´ı tˇri hodnoty xi a yi z druh´e ˇc´asti tabulky 7.3.1, tj. x0 = −2, x1 = −1, x2 = 0 a y0 = −1, y1 = 1.2508, y2 = 1.3112. Protoˇze f (x, y) = x2 − 0.2y, vypoˇc´ıt´ame f0 = (−2)2 − 0.2 · (−1) = 4.2000, f1 = (−1)2 − 0.2 · 1.2508 = 0.7498 a f2 = (0)2 − 0.2 · 1.3112 = −0.2622. Pomoc´ı vzorce (7.4.5) pak provedeme vˇsechny dalˇs´ı v´ ypoˇcty: y3 = −1.3112 +
h (23f2 12
− 16f1 + 5f0 ) = 1.5588,
f3 = 12 − 0.2 · 1.5588 = 0.6882, y4 = 1.5588 +
h (23f3 12
− 16f2 + 5f1 ) = 3.5400,
f4 = 22 − 0.2 · 3.5400 = 3.2920, y5 = 3.5400 +
h (23f4 12
− 16f3 + 5f2 ) = 8.8227.
V´ysledky jsou shrnuty v tabulce 7.4.1. Adamsovy-Moultonovvy vzorce Adamsovy-Moultonovy vzorce jsou implicitn´ı v´ıcekrokov´e metody, kter´e lze odvodit podobn´ ym postupem jako Adamsovy-Bashforthovy vzorce, tj. provede se integrace interpolaˇcn´ıho polynomu sestaven´eho pro derivaci ˇreˇsen´ı diferenci´aln´ı rovnice. Nyn´ı se integruje na intervalu xi−1 , xi a pak se posune indexov´an´ı. Uved’me pouze pˇrehled vzorc˚ u: yi+1 = yi + hfi+1 ,
(7.4.7)
- 146 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
h yi+1 = yi + (fi+1 + fi ), 2 h yi+1 = yi + (5fi+1 + 8fi − fi−1 ), 12 h yi+1 = yi + (9fi+1 + 19fi − 5fi−1 + fi−2 ). 24
(7.4.8) (7.4.9) (7.4.10)
Napˇr´ıklad (7.4.9) pˇredsatvuje dvoukrokovou, jednobodouvou metodu tˇret´ıho ˇra´du. Pouˇzit´ı uveden´ych vzorc˚ u vyˇzaduje v kaˇzd´em kroku ˇreˇsit rovnici pro nezn´amou usob ˇreˇsen´ı, kter´ y jsme zm´ınili v u ´ vodu, rozebereme aˇz v dalˇs´ım yi+1 . Iteraˇcn´ı zp˚ odstavci. Zde si ukaˇzme jednoduˇzˇs´ı postup, kdy se do zvolen´eho vzorce dosad´ı prav´a strana diferenci´aln´ı rovnice f a pak se vyj´adˇren´ım yi+1 vzorec pˇrevede na explicitn´ı tvar. Takto m˚ uˇzeme ovˇsem postupovat pouze v pˇr´ıpadech, kdy prav´a strana f nen´ı pˇr´ıliˇs sloˇzit´a funkce”. ” Tabulka 7.4.1: Adamsova-Bashforthova a Adamsova-Moultonova metoda. Ad.-Bash. 3. ˇra´du i
xi
yi
0 −2 1 −1 2 0 3 1 4 2 5 3
−1 1.2508 1.3112 1.5588 3.5400 8.8227
Ad.-Moul. 3. ˇra´du
|yi − y(xi )| 0 0.0001 0.0001 0.1679 0.2411 0.3060
i
xi
yi
0 −2 1 −1 2 0 3 1 4 2 5 3
−1 1.2508 1.2929 1.3613 3.2628 8.4773
|yi − y(xi )| 0 0.0001 0.0183 0.0296 0.0362 0.0394
Pˇ r´ıklad 7.4.2. Poˇc´ateˇcn´ı u ´ lohu (7.1.4) ˇreˇste na intervalu −2, 3 pomoc´ı dvoukrokov´e Adamsova-Moultonovy metody pro h = 1. V´ ysledky porovnejte s pˇresn´ym ˇreˇsen´ım. ˇ sen´ı: Do prav´e strany vzorce (7.4.9) dosad´ıme fi+1 = f (xi+1 , yi+1) = x2i+1 − Reˇ 0.2yi+1 : yi+1 = yi +
h (5(x2i+1 − 0.2yi+1) + 8fi − fi−1 ). 12
- 147 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
Vyjadˇr´ıme yi+1: yi+1 = 12(yi +
h (5x2i+1 + 8fi − fi−1 ))/(12 + h). 12
Tento pˇrepis Adamsova-Moultonova vzorce pouˇzijeme ve v´ypoˇctech. Jako poˇc´ateˇcn´ı u ´ sek vezmeme hodnotu vypoˇc´ıtanou metodou RK4, tj. pouˇzijeme xi a yi z prvn´ıch dvou ˇra´dk˚ u tabulky 7.3.1. V´ ysledky jsou shrnuty v tabulce 7.4.1. Algoritmus prediktor-korektor Princip odvozen´ı algoritmu uk´aˇzeme na dvoukrokov´em Adamsovˇe-Moultonovˇe vzorci: yi+1 = yi +
h (5f (xi+1 , yi+1) + 8fi − fi−1 ). 12
Pˇredpokl´adejme, ˇze zn´ame hodnoty yi, yi−1 , yi−2 a chceme vypoˇc´ıtat yi+1 . Na uveden´y vzorec m˚ uˇzeme pohl´ıˇzet jako na rovnici v iteraˇcn´ım tvaru (2.4.1), pro jej´ıˇz vyˇreˇsen´ı je pˇrirozen´e pouˇz´ıt metodu prost´ ych iterac´ı (2.4.3). Tato u ´vaha vede na n´asleduj´ıc´ı v´ ypoˇcetn´ı postup: Pro i = 1, . . . , n − 1 vypoˇcti yi+1 takto: 0 , (a) urˇci poˇca´teˇcn´ı odhad yi+1
(b) poˇc´ateˇcn´ı odhad zpˇresˇ nuj pomoc´ı iterac´ı: k+1 = yi + yi+1
h k (5f (xi+1 , yi+1 ) 12
+ 8fi − fi−1 ),
k = 0, 1, 2, . . . .
Krok (a) se naz´ yv´a prediktor, krok (b) je korektor. Varianta algoritmu prediktorkorektor uveden´a n´ıˇze pouˇz´ıv´a explicitn´ı Adams˚ uv-Bashforth˚ uv vzorec tˇret´ıho ˇra´du (7.4.5) jako prediktor. V korektoru se pak vykon´av´a jedin´e iteraˇcn´ı zpˇresnˇen´ı. S ohledem na struˇcn´y z´apisu algoritmu pouˇzijeme symbol pˇriˇrazen´ı :=” pro ” zmˇenu hodnoty promˇenn´e, coˇz umoˇzn ˇ uje vynechat iteraˇcn´ı index k. Dost´av´ame n´asleduj´ıc´ı algoritmus:
- 148 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
y0 = c; y1 , y2 vypoˇcti pomoc´ı jednokrokov´e metody; Pro i = 2, . . . , n − 1 vypoˇcti yi+1 takto: (P ) yi+1 := yi +
h (23fi 12
− 16fi−1 + 5fi−2 ),
(E) fi+1 := f (xi+1 , yi+1), (C) yi+1 := yi +
h (5fi+1 12
+ 8fi − fi−1 ),
(E) fi+1 := f (xi+1 , yi+1).
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
(7.4.11)
Uveden´y algoritmus je metoda tˇr´ıkrokov´ a, dvoubodov´a a tˇret´ıho ˇra´du, pro pro niˇz se pouˇz´ıv´a oznaˇcen´ı PECE. Podobnˇe se pouˇz´ıv´a oznaˇcen´ı PE(CE)k , P(EC)k nebo P(EC)k E, pro varianty algoritmu prediktor-korektor s k vnitˇrn´ımi kroky a s r˚ uznou organizac´ı dalˇs´ıch v´ ypoˇct˚ u. Pˇ r´ıklad 7.4.3. Poˇc´ateˇcn´ı u ´ lohu (7.1.4) ˇreˇste na intervalu −2, 3 pomoc´ı PECE metody (7.4.11) s krokem h = 1. V´ ysledky porovnejte s pˇresn´ym ˇreˇsen´ım. ˇ sen´ı: Poˇc´ateˇcn´ı u Reˇ ´sek vypoˇc´ıt´ame metodou RK4. Jsou to prvn´ı tˇri hodnoty xi a yi z druh´e ˇc´asti tabulky 7.3.1, tj. x0 = −2, x1 = −1, x2 = 0 a y0 = −1, y1 = 1.2508, y2 = 1.3112. Pomoc´ı tˇechto hodnot vypoˇc´ıt´ame f0 = 4.2000, f1 = 0.7498 ypoˇct˚ u podle (7.4.11) uv´ad´ıme pouze poˇrad´ı v jak´em a f2 = −0.2622. U dalˇs´ıch v´ vznikaj´ı jednotliv´e hodnoty: y3 := 1.5588, f 3 := 0.6882, y3 := 1.3607, f 3 := 0.7279, y4 := 3.4178, f 4 := 3.3164, y4 := 3.2496, f 4 := 3.3501, y5 := 8.5908, f 5 := 7.2818, y5 := 8.4564, f 5 := 7.3087. V´ysledky jsou v tabulce 7.4.2.
- 149 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
Tabulka 7.4.2: Algoritmus prediktor-korektor. i
xi
yi
0 −2 1 −1 2 0 3 1 4 2 5 3
−1 1.2508 1.3112 1.3607 3.2496 8.4564
|yi − y(xi )| 0.0000 0.0001 0.0001 0.0302 0.0493 0.0603
Pozn´ amka Z porovn´ an´ım v´ysledk˚ u v tabulk´ ach 7.4.1 a 7.4.2 je vidˇet, ˇze implicitn´ı metoda je podstatnˇe pˇresnˇejˇs´ı neˇz explicitn´ı metoda stejn´eho ˇra´du. Algoritmus prediktorkorektor pˇredstavuje nepˇresnou” implicitn´ı metodu, takˇze tak´e dosaˇzen´ a pˇresnost ” je o nˇeco menˇs´ı.
Kontroln´ı ot´ azky Ot´azka 1. V ˇcem spoˇc´ıv´a hlavn´ı pˇr´ınos v´ıcekrokov´ ych metod oproti metod´am jednokrokov´ ym? Ot´azka 2. Jak se odvozuj´ı Adamsovy-Bashforthovy a Adamsovy-Moultonovy vzorce? Ot´azka 3. Na ˇcem je zaloˇzena konstrukce algoritmu prediktor-korektor?
´ Ulohy k samostatn´ emu ˇreˇsen´ı 1. Diferenci´aln´ı rovnici y = y(1 + sin x) − x3 , y(1) = 0 ˇreˇste na intervalu 1, 2 s krokem h = 0.2 pomoc´ı Adamsovy-Bashforthovy metody ˇctvrt´eho ˇra´du. 2. Vzorec Adamsovy-Moultonovvy metody ˇctvrt´eho ˇra´du pˇrepiˇste na explicitn´ı tvar pro diferenci´aln´ı rovnici z pˇredchoz´ı u ´ lohy a proved’te v´ ypoˇcet. 3. Sestavte algoritmy prediktor korektor PECE a PECEC zaloˇzen´e na Adamsov´ ych-Bashforthov´ ych a Adamsov´ych-Moultonov´ych vzorc´ıch ˇctvrt´eho ˇra´du a vyˇreˇste diferenci´aln´ı rovnici z prvn´ı u ´ lohy.
- 150 -
ˇ ATE ´ CN ˇ ´I ULOHY ´ 7. ODR – POC
Numerick´e metody
V´ ysledky u ´loh k samostatn´ emu ˇreˇsen´ı 1. Poˇc´ateˇcn´ı u ´sek vypoˇc´ıt´ame metodou RK4, pak y4 = −4.6497, y5 = −8.5164. h (−9x3i+1 + 19fi −5fi−1 + fi−2 ))/(8 −3h(1 + sin xi+1 )). 2. Odvod´ıme yi+1 = 8(yi + 24
Pomoc´ı y0 , y1, y2 z RK4 vypoˇc´ıt´ame: y3 = −2.3270, y4 = −4.6615, y5 = −8.5396. 3. Algoritmus PECEC m´a n´asleduj´ıc´ı podobu: y0 = c; y1 , y2 , y3 vypoˇcti pomoc´ı jednokrokov´e metody (RK4); Pro i = 3, . . . , n − 1 vypoˇcti yi+1 takto: (P ) yi+1 := yi +
h (55fi 24
− 59fi−1 + 37fi−2 − 9fi−3 ),
(E) fi+1 := f (xi+1 , yi+1), (C) yi+1 := yi +
h (9fi+1 24
+ 19fi − 5fi−1 + fi−2 ),
(E) fi+1 := f (xi+1 , yi+1), (C) yi+1 := yi +
h (9fi+1 24
+ 19fi − 5fi−1 + fi−2 ).
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
Algoritmus PECE dostaneme vynech´an´ım posledn´ıho kroku (C). Poˇc´ateˇcn´ı u ´ sek vypoˇc´ıt´ame opˇet metodou RK4. V´ ysledky podle PECE: y4 = −4.6581, y5 = −8.5342; v´ ysledky podle PECEC: y4 = −4.6594, y5 = −8.5360.
Shrnut´ı lekce Uk´azali jsme princip numerick´eho ˇreˇsen´ı poˇc´ateˇcn´ıch u ´ loh pro obyˇcejn´e diferenci´aln´ı rovnice (prvn´ıho ˇra´du). Metody tohoto typu lze nal´ezt ve standardn´ıch knihovn´ach poˇc´ıtaˇcov´ ych program˚ u. Popisovan´a problematika je pomˇernˇe obs´ahl´a a zahrnuje napˇr´ıklad ˇr´ızen´ı dosaˇzen´e pˇresnosti nebo speci´aln´ı metody pro diferenci´aln´ı rovnice, jejichˇz ˇreˇsen´ı je citliv´e na zaokrouhlovac´ı chyby. ´ Ulohy zapsan´e pomoc´ı diferenci´aln´ıch rovnic (r˚ uzn´ych typ˚ u nejen ODR) se vyskytuj´ı velice ˇcasto pˇri modelov´an´ı fyzik´aln´ıch nebo technick´ych u ´ loh. Jejich ˇreˇsen´ı na poˇc´ıtaˇci je atraktivn´ı oblast´ı v´ yzkumu, kter´ y v posledn´ıch desetilet´ıch proch´az´ı bouˇrliv´ ym v´ yvojem.
- 151 -
