Pokročilé numerické metody II
1. přednáška
Základní pojmy
Jednokrokové metody
Eulerovy metody
Jiří Zelinka
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 1 / 15
Literatura
Rektorys, K.: Variační metody v inženýrských problémech
a v problémech matematické fyziky. Vyd. 6., opr. české 2.
Praha: Academia, 1999. 602 s. ISBN 8020007148.
Vitásek, E.: Základy teorie numerických metod pro řešení
diferenciálních rovnic. 1. vyd. Praha: Academia, 1994. 409
s. ISBN 8020002812.
Ralston, A.: Základy numerické matematiky. Translated
by Milan Práger - Emil Vitásek. 2. čes. vyd. Praha:
Academia, 1978.
Mathews, J.H., Fink, K.D.: Numerical methods using
MATLAB, Pearson Prentice Hall, 2003
Stoer, J., Bulirsch R.: Introduction to Numerical Analysis,
Spriger, 1992
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 2 / 15
Obsah přednášky
1 Obyčejné diferenciální rovnice
jednokrokové metody
vícekrokové metody
okrajové úlohy
– metoda střelby
– diferenční metody
2 Variační metody (ODR i PDR)
3 Diferenční metody pro PDR
eliptické rovnice
časově závislé rovnice
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 3 / 15
Následující část je převzata z 6. kapitoly výukových materiálů
Numerické metody Jaromíra Kuben a Pavlíny Račkové,
Univerzita obrany.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 4 / 15
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
568
Kapitola 6
Numerické řešení obyčejných diferenciálních
rovnic
Cíle
Po prostudování této kapitoly budete schopni vysvětlit:
co je to obyčejná diferenciální rovnice (ODR) prvního řádu v implicitním
resp. explicitním tvaru a co je její řešení,
jak vypadá Cauchyova počáteční úloha pro ODR prvního řádu,
proč je nutné hledat řešení počáteční úlohy numericky a jaký je princip
numerického řešení,
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 569
jaký je rozdíl mezi jednokrokovými a vícekrokovými metodami,
jaký je obecný tvar jednokrokové explicitní a implicitní metody,
co je to lokální a globální chyba, řád, konvergence a A-stabilita jednokrokové
metody,
co jsou to metody Rungeho-Kutty a jaký je jejich princip,
jaký je obecný tvar lineární vícekrokové explicitní a implicitní metody,
co je to lokální a globální chyba, řád, konvergence a A-stabilita lineární
vícekrokové metody,
co jsou Adamsovy-Bashforthovy a Adamsovy-Moultonovy metody a metody
zpětného derivování a jaký je princip jejich odvození,
co jsou metody prediktor-korektor,
čemu se říká tuhé problémy a na co je nutné dbát při jejich řešení,
jak se výsledky pro jednu ODR prvního řádu přenesou na soustavy ODR
prvního řádu,
jak přepíšeme jednu ODR vyššího řádu na soustavu ODR prvního řádu.
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 570
Jde o jednu z nejdůležitějších úloh numerické matematiky, protože diferenciální rovnice
patří k nejvýznamnějším matematickým modelům reálných problémů. Přitom pouze
malá část obyčejných diferenciálních rovnic má řešení, které lze vyjádřit v uzavřeném
tvaru, tj. pomocí elementárních funkcí.
6.1 Obyčejné diferenciální rovnice 1. řádu
6.1.1 Základní vlastnosti obyčejných diferenciálních rovnic 1. řádu
Připomeneme si základní poznatky z tzv. kvalitativní teorie obyčejných diferenciálních
rovnic. Podrobněji viz [34].
Rovnice F.x; y; y0
/ D 0 se nazývá obyčejná diferenciální rovnice prvního řádu v implicitním
tvaru s neznámou funkcí y.x/ nezávisle proměnné x. Příkladem takové
rovnice je
ln.x2
yy0
C y2
/ .y0
/3
C x sin.xyy0
/
p
y y0 D 0:
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 571
Pro vyšetřování rovnic je důležitý případ, kdy se dá osamostatnit první derivace y0
, tj.
rovnice má tvar
y0
D f .x; y/: (6.1)
O takové rovnici říkáme, že je v explicitním tvaru. O funkci f .x; y/ předpokládáme,
že je definována na nějaké otevřené množině D R2
.
Příklady takových rovnic jsou
y0
D x2
C y2
; y0
D
y
x
; y0
D y2
; y0
D sin xy; y0
D
p
x2 y2 apod.
Řešením je funkce y.x/ definovaná na intervalu I, která splňuje rovnici (6.1). Musí tedy
platit, že Œx; y.x/ 2 D pro x 2 I (jinými slovy, graf funkce y.x/ leží v množině D —
viz obr. 6.1 a)) a y0
.x/ D f .x; y.x// pro každé x 2 I.
V každém bodě Œx; y.x/ grafu řešení tedy platí, že tečna ke grafu v tomto bodě má
směrnici rovnou číslu f .x; y.x//. Protože funkci f .x; y/ známe, můžeme sestrojit
(teoreticky) v každém bodě Œx; y vázaný vektor, jehož směrnice je f .x; y/. Graf řešení
se musí v každém bodě, kterým prochází, dotýkat tohoto vektoru. Tyto vektory tvoří
směrové pole diferenciální rovnice. Na obr. 6.2, 6.3 a 6.4 jsou znázorněna směrová pole
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 572
x
y
0
I
y.x/
D
a) Řešení diferenciální rovnice
x
y
y0
0 x0
y.x/
b) Cauchyova počáteční úloha
Obr. 6.1
a několik řešení šesti diferenciálních rovnic. Všimněte si, že u rovnic na obr. 6.3 a)
až 6.4 b) není množina D rovna celé rovině R2
(určete definiční obory pravých stran
těchto rovnic).
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 573
a) y0
D x2
C y2
b) y0
D sin.xy/
Obr. 6.2: Směrová pole a řešení diferenciálních rovnic — část 1
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 574
a) y0
D 2e x
.1 C 1=x/y b) y0
D 1=.x C 2y/
Obr. 6.3: Směrová pole a řešení diferenciálních rovnic — část 2
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 575
a) y0
D
p
9 x2 y2 b) y0
D
p
y x2=2 C 3
Obr. 6.4: Směrová pole a řešení diferenciálních rovnic — část 3
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 576
Rovnice y0
D f .x; y/ má obvykle nekonečně mnoho řešení. Hledáme tedy řešení,
které splňuje nějakou dodatečnou podmínku. Tato podmínka by měla jednoznačně určit
jediné řešení.
My budeme hledat řešení, jehož graf prochází zadaným bodem Œx0; y0. Chceme tedy,
aby platila tzv. počáteční podmínka y.x0/ D y0 — viz obr. 6.1 b). Úloha, kterou
budeme řešit, má tudíž tvar
y0
.x/ D f .x; y.x//;
y.x0/ D y0:
(6.2)
Říká se jí obvykle Cauchyova počáteční úloha.
Připomeňme si, jaké vlastnosti funkce f .x; y/ zaručují, že Cauchyova počáteční
úloha (6.2) má řešení a že toto řešení je jediné.
1) Z teorie obyčejných diferenciálních rovnic je známo, že spojitost funkce f .x; y/
zaručuje existenci řešení počáteční úlohy.
Toto řešení ale nemusí být jediné. Pak říkáme, že v bodě Œx0; y0 je porušena jednoznačnost.
Protože všechna řešení počáteční úlohy (6.2) musí mít v bodě Œx0; y0
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 577
společnou tečnu t, musí dojít k jakémusi „rozvětvení“ — viz obr. 6.5. Numerické
hledání řešení je pak velmi problematické.
x
y
x00
y0
t
Obr. 6.5: Porušení jednoznačnosti v bodě Œx0; y0
Např. počáteční úloha y0
D 2
p
jyj, y.0/ D 0, má spojitou pravou stranu f .x; y/ D
D 2
p
jyj v množině D D R2
. Uvedená úloha má ale nekonečně mnoho řešení.
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 578
Jedním je nulová funkce y.x/ Á 0, dalšími řešeními jsou funkce
yc D
´
.x c/2
; x = c;
0; x < c;
kde c 2 RC
0 je libovolné číslo. Řešení jsou znázorněna na obr. 6.6. Že jde opravdu
o řešení, se snadno ověří dosazením do rovnice.
2) Standardní vlastností, která zaručuje, že počáteční úloha (6.2) má jediné řešení,
je Lipschitzova podmínka (existuje konstanta L > 0 taková, že je v množině D
splněna nerovnost jf .x; y1/ f .x; y2/j 5 Ljy1 y2j; stačí, aby tato vlastnost
platila lokálně). K platnosti této podmínky stačí, aby v okolí bodu Œx0; y0 byla spojitá
parciální derivace @
@y
f .x; y/. Tuto podmínku rovnice, se kterými se v aplikacích
setkáváme, většinou splňují.
Titulní strana Copyright Předmluva Obsah Kap1 Kap2 Kap3 Kap4 Kap5 Kap6 Kap7 Testy Literatura Rejstřík Tiráž
Obsah Jdi na stranu J J I I J I Celá obrazovka/Okno Zavřít
Numerické řešení obyčejných diferenciálních rovnic 579
x
y
1 2;5 4 c0
x2
.x 1/2
.x 2;5/2
.x 4/2
.x c/2
Obr. 6.6: Řešení úlohy y0
D 2
p
jyj, y.0/ D 0
Systémy diferenciálních rovnic
y′
1 = f1(x, y1, . . . , yn), y1(x0) = y10
y′
2 = f2(x, y1, . . . , yn), y2(x0) = y20
...
y′
n = fn(x, y1, . . . , yn), yn(x0) = yn0
Zápis jednou rovnicí
Y ′
= F(x, Y ), Y (x0) = Y0,
Y , F jsou vektory.
Rovnice vyšších řádů
y(n)
= f (x, y, y′
, y′′
, . . . , y(n−1)
),
jsou dány počáteční podmínky y(x0), y(x0),. . . , y(n−1)
(x0).
Lze převést na systém DR.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 5 / 15
Numerické řešení počátečního problému
y′
= f (x, y), y(x0) = y0
Numerické řešení hledáme v bodech (uzlech)
x0 < x1 < · · · < xn, hi = xi+1 − xi ,
zpravidla používáme ekvidistantní uzly, tj. hi = h, ∀i.
Označení: yi je hodnota numerického řešení v bodě xi , tj.
yi ≈ y(xi ).
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 6 / 15
Jednokrokové metody
Obecná explicitní jednokroková metoda
yi+1 = yi + h · Φ(xi , h, yi )
Funkce Φ se nazývá přírůstková funkce, závisí i na f .
Poznámka: Φ může záviset i na yi+1, pak se jedná implicitní
metodu.
Lokální diskretizační chyba
Jedná se o chybu metody v jednom kroku, pokud bychom
použili přesné hodnoty řešení. Označení: ltei
ltei = y(xi+1) − y(xi ) − hΦ(xi , y(xi ), h, y(xi+1))
Lokální chyba
Bodem (xi , yi ) prochází nějaké řešení ui diferenciální rovnice,
lokální chyba lei je dána chybou metody pro ui :
lei = ui (xi+1) − yi+1
Globální chyba ei = y(xi ) − yi .
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 7 / 15
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 8 / 15
Řád a konvergence metody
Číslo p nazýváme řádem metody, jestliže ltei je O(hp+1
).
Věta
Jestliže ltei = O(hp+1
) a funkce Φ je spojitá a lipschitzovská
v y, tj. existuje konstanta M, že
|Φ(x, h, ˜y) − Φ(x, h, ˆy)| ≤ M|˜y − ˆy|,
pak ei = O(hp
).
Lemma: Jestliže pro posloupnost a0, a1, . . . platí
|ai+1| ≤ (1 + λ)|ai | + B, λ, B > 0, pak |an| ≤ eλn|a0| + B eλn−1
λ .
Metoda se nazývá konvergentní v bodě x(= xn), jestliže pro
h → 0, n → ∞, nh = xn − x0 platí yn → y(x).
Metoda se nazývá konvergentní na intervalu [a, b], (zpravidla
a = x0), jestliže je konvergentní ∀x ∈ [a, b].
Poznámka: konvergence zpravidla plyne z řádu globální chyby.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 9 / 15
Explicitní Eulerova metoda
Taylorův rozvoj:
y(x +h) = y(x)+hy′
(x)+O(h2
) = y(x)+hf (x, y(x))+O(h2
),
tedy
y(xi+1) = y(xi ) + hf (xi , y(xi )) + O(h2
).
Zanedbáním chybového členu a použitím přibližných hodnot
dostaneme Eulerovu explicitní metodu:
yi+1 = yi + hf (xi , yi )
Současně vidíme, že ltei = O(h2
), takže metoda je řádu 1.
Dá se ukázat, že i lei = O(h2
) a ltei − lei = O(h3
).
Dále je Φ(x, h, y) = f (x, y), takže pokud funkce f splňuje
předpoklady Picardovy věty o existenci a jednoznačnosti řešení
rovnice, pak globální chyba je O(h) a metoda je konvergentní.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 10 / 15
Implicitní Eulerova metoda
Taylorův rozvoj:
y(x) = y(x + h) − hy′
(x + h) + O(h2
) =
= y(x + h) − hf (x + h, y(x + h)) + O(h2
),
tedy
y(xi ) = y(xi+1) − hf (xi+1, y(xi+1)) + O(h2
).
Zanedbáním chybového členu, převedením na druhou stranu
rovnice a použitím přibližných hodnot dostaneme Eulerovu
implicitní metodu:
yi+1 = yi + hf (xi , yi+1)
Opět platí ltei = O(h2
), podobně jsou další vlastnosti (řád
globální chyby a konvergence metody) stejné jako u explicitní
metody.
V každém kroku Eulerovy implicitní metody ovšem musíme
řešit nelineární rovnici.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 11 / 15
Lichoběžníková metoda
Integrální tvar ODR:
y(x) = y0 +
x
x0
f (t, y(t))dt
y(xi+1) = y(xi ) +
xi+1
xi
f (t, y(t))dt
Lichoběžníkové pravidlo pro numerický výpočet integrálu:
yi+1 = yi +
1
2
h[f (xi , yi ) + f (xi+1, yi+1)]
Chyba lichoběžníkového pravidla je O(h3
), takže i lokální
diskretizační chyba pro lichoběžníkovou metodu je O(h3
), tedy
metoda je řádu 2 a je konvergentní.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 12 / 15
Metoda prediktor–korektor
Iterační řešení implicitní metody
Prediktor – explicitní Eulerova metoda – počáteční iterace:
y
(0)
i+1 = yi + hf (xi , yi )
Korektor – implicitní Eulerova metoda – zpřesňování iterace:
y
(j+1)
i+1 = yi + hf (xi+1, y
(j)
i+1)
– podobně s lichoběžníkovou metodou
Metoda konverguje pro dostatečně malé h.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 13 / 15
Stabilita
Zkoumáme, co se děje pro konstantní h a n → ∞.
Pro testovací úlohu y′
= λy, y(0) = 1, λ ∈ C, R(λ) < 0 je
řešení y(x) = eλx
. Pro něj platí |y(xn+1)| < |y(xn)|, proto
požadujeme aby
|yn+1| < |yn|.
Pro obecnou jednokrokovou metodu a testovací úlohu
dostaneme vztah
yn+1 = R(z)yn,
kde R se nazývá funkce stability, z = λh. Zajímá nás tedy
podmínka |R(z)| < 1.
Oblast (absolutní) stability: Ω = {z ∈ C : |R(z)| < 1}.
Interval stability: reálná část oblasti stability.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 14 / 15
Stabilita
Metoda se nazývá absolutně stabilní (A–stabilní), pokud Ω
obsahuje všechna záporná reálná čísla.
Metoda se nazývá L–stabilní, jestliže je A–stabilní a
lim
|z|→∞
R(z) = 0.
Oblast stability pro explicitní Eulerovu metodu: |z − (−1)| < 1.
Implicitní Eulerova metoda a lichoběžníková metoda jsou
absolutně stabilní, implicitní Eulerova metoda je dokonce
L–stabilní.
Jiří Zelinka Pokročilé numerické metody II, 1. přednáška 15 / 15