Numerické metody pro řešení počátečních
úloh pro ODR – jednokrokové metody
Mnoho matematických modelů systémů, jejichž stavové proměnné se podle jistého zákona
mění v závislosti na okamžitém stavu systému, lze popsat pomocí obyčejné diferenciální
rovnice, resp. soustav diferenciálních rovnic, prvního řádu. Řešení, pokud existuje,
není jediné. Abychom zaručili jednoznačnost řešení, požadujeme ještě splnění tzv. počáteční
podmínky.
Formulace: Hledáme řešení y = y(x) rovnice (1) s počáteční podmínkou (2)
y (x) = f(x, y(x))(1)
y(x0) = y0.(2)
Smysl: Analyticky lze spočítat jen velmi malou skupinu počátečních úloh pro ODR.
Proto je tak důležité numerické řešení.
Princip: Základem metod je diskretizace proměnných. Přibližné řešení se nekonstruuje
jako spojitá funkce, ale nagenerujeme body x0, x1, x2, . . . a určujeme čísla y0, y1, y2, . . .,
která aproximují y(x0), y(x1), y(x2), . . ..
Poznámka: Body sítě x0, x1, x2, . . . nemusí být ekvidistantní:
xi+1 = xi + hi.
Platí-li: hi = h ∀i mluvíme o metodě s konstantním krokem
(ekvidistantní síť)
Neplatí-li: hi = h ∀i mluvíme o metodě s proměnným krokem
Poznámka: Aproximace yn hodnoty přesného řešení y(xn) v bodě xn se počítá z hodnot
přibližného řešení v předchozích uzlech.
Počítáme-li yn+1 pouze pomocí hodnoty yn mluvíme o jednokrokové metodě.
Počítáme-li yn+1 pomocí více předchozích hodnot yn, yn−1, . . . mluvíme
o vícekrokové metodě.
1
Jednokrokové metody
Nejjednodušší metodou je Eulerova metoda.
Princip:
y0 . . . je dáno (počáteční podmínka)
y1 . . . počítáme extrapolací z hodnoty y0, přičemž se na intervalu
x0, x1 řešení aproximuje přímkou, která prochází bodem [x0, y0]
a má směrnici y = f(x0, y0).
Ta má rovnici y = y0 + (x − x0)f(x0, y0).
Tj. pro x1 dostáváme:
y1 = y0 + (x1 − x0)
h0
·f(x0, y0).
Obecně dostaneme rekurentní vztah:
yn+1 = yn + hn · f(xn, yn), n = 0, 1, 2, . . .
Geometricky:
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
6
x
y přesné řešení . . . y(x)
x0 x1 x2 x3 x4
y0
y1
y2
y3
y4
y(x1)
y(x2)
y(x3)
y(x4)
2
Poznámky:
1. Eulerovu metodu můžeme chápat také tak, že hodnotu y(xn+1) = y(xn + hn) aproximujeme
pomocí Taylorova polynomu 1 stupně pro funkci y v bodě xn:
y(xn+1) ≈ y(xn) + hny (xn) = y(xn) + hnf(xn, y(xn)).
2. Také ji lze chápat tak, že diferenciální rovnici y = f(x, y) nahradíme diferenční
rovnicí
yn+1 − yn
hn
= f(xn, yn) n = 0, 1, 2, . . .
3
Příklad: Řešte úlohu
y = x − y,
y(0) = 1
na intervalu 0; 0,6 s konstantními
kroky h = 0, 2 a h = 0, 1.
(Přesné řešení: y(x) = 2e−x
+ x − 1).
Řešení: Použijeme rekurentní vztah:
yn+1 = yn + h · f(xn, yn).
h = 0, 2 h = 0, 1
xn
přesné
y(xn) yn en yn en
0 1,000 1,000 0,000 1,000 0,000
0,1 0,910 0,900 0,010
0,2 0,837 0,800 0,037 0,820 0,017
0,3 0,782 0,758 0,024
0,4 0,741 0,680 0,061 0,712 0,029
0,5 0,713 0,681 0,032
0,6 0,698 0,624 0,074 0,663 0,035
−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0
0.5
1
1.5
x
y
. . . přesné řešení
. . . řešení pro h = 0, 1
. . . řešení pro h = 0, 2
Poznámka: 1) Vidíme, že je chyba úměrná h,
2) Chyba s rostoucím x vzrůstá.
4
Obecná jednokroková metoda
Eulerova metoda je sice velmi jednoduchá, ale k dosažení určité přesnosti musíme používat
velmi malé kroky hi. Chceme-li jednokrokovou metodu vyššího řádu, musíme se zříci
linearity, tj.
yn+1 = yn + \hn · Φ(xn, yn, hn, f) n = 0, 1, 2, . . .
Metody Taylorova typu:
Hodnotu y(xn+1) budeme aproximovat pomocí Taylorova rozvoje vyššího řádu
(1. řádu = Eulerova metoda), tj.
y(xn+1) = y(xn + hn) = y(xn) + hny (xn) +
h2
n
2!
y (xn) + . . . +
hp
n
p!
y(p)
(xn)(3)
Derivace y v bodě xn lze určit postupným derivováním funkce f.
y = f(x, y(x))
y = fx + fy · y
=f(x,y(x))
=
∂f
∂x
fx
+
∂f
∂y
fy
.
dy
dx
y
Obecně lze odvodit rekurenci:
y = f(x, y(x))
y(r+1)
= f(r)
(x, y(x)) = f(r−1)
x (x, y(x)) + f(r−1)
y (x, y(x)) · f(x, y(x)) r = 1, 2, . . .
(4)
Zbývá jen dosadit (4) za derivace v (3).
5
Příklad: Odvoďte metodu Taylorova typu 2.řádu pro řešení úlohy:
y = x − y,
y(0) = 1
na intervalu 0; 0, 6 s konstantním krokem h = 0, 2.
(Přesné řešení: y(x) = 2e−x
+ x − 1).
Řešení:
f(x, y) = x − y
f (x, y) = fx + fy · f = 1 + (−1) · f(x, y) = 1 − x + y.
Dostáváme rekurentní vztah:
yn+1 = yn + hn(xn − yn) +
1
2
h2
n(1 − xn + yn)
xn
přesné
y(xn) yn h(xn − yn)
h2
2
(1 − xn + yn) en
0 1,000 1,000 -0,200 0,040 0,000
0,2 0,837 0,840 -0,128 0,033 -0,003
0,4 0,741 0,745 -0,069 0,027 -0,004
0,6 0,698 0,703 -0,005
Poznámka: Vidíme, že metoda Taylorova typu 2. řádu pro h = 0, 2 dává přesnější
výsledky než Eulerova metoda s h = 0, 1.
6
Metody Runge-Kuttova typu
• Univerzálnější metody než metody Taylorova typu.
• Vychází také z Taylorova polynomu, ale nepoužívá se ho přímo, aby nebylo nutné
explicitně vyjadřovat derivace funkce f = f(x, y(x)) a počítat jejich hodnoty. Hledaná
aproximace je kombinací několika hodnot funkce f vypočítaných v několika
strategicky volených bodech (x, y) na intervalu xn, xn+1 .
Poznámka: Těchto metod je velké množství!
Ukážeme si odvození dvou metod tohoto typu s geometrickou interpretací.
Použijeme následující úvahy:
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
6
x
y
x1x0 x0 + x1
2
M0
P
M1
Věta:
Nechť oblouk M0 M1 je částí paraboly. Potom platí:
1. Tečna v bodě P je rovnoběžná s tětivou M0 M1.
2. Směrnice tětivy M0 M1 je aritmetickým průměrem směrnic tečen v M1 a M2.
7
Důkaz: Rovnice paraboly (polynomu 2.stupně): y − b = c(x − a)2
y = c(x − a)2
+ b
⇒
y = 2c(x − a)
0 2 4 6 8
−1
0
1
2
3
4
5
6
x
y
a
b
y − b = c(x − a)2
1. Směrnice tečny v bodě P:
y (
x0 + x1
2
) = 2c(
x0 + x1
2
− a) = c(x0 + x1 − 2a)
Směrnice tětivy M0M1 je:
y(x1) − y(x0)
x1 − x0
=
c(x1 − a)2
+ b − c(x0 − a)2
− b
x1 − x0
=
=
cx2
1 − 2acx1 + a2
c + b − cx2
0 + 2acx0 − a2
c − b
x1 − x0
=
= c
x2
1 − x2
0 − 2a(x1 − x0)
x1 − x0
= c(x1 + x0 − 2a).
2. Směrnice tečny v bodě M0 je:
y (x0) = 2c(x0 − a)
Směrnice tečny v bodě M1 je:
y (x1) = 2c(x1 − a)
Jejich aritmetický průměr:
y (x0) + y (x1)
2
=
2c(x0 − a) + 2c(x1 − a)
2
=
= c(x0 − a + x1 − a) = c(x0 + x1 − 2a).
8
Nyní použijeme vlastnost 1)
Známe souřadnice bodu M0. Jestliže bychom znali y-souřadnici bodu P, pak stačí
udělat tečnu a bodem M0 vést rovnoběžku a dostaneme y-souřadnici bodu M1. My
ale y-souřadnici bodu P neznáme (obecně funkce y = y(x) nemusí být parabola, to
je jen naše aproximace), takže ji vyjádříme přibližně. Bod P nahradíme bodem P ,
který má stejnou x-ovou souřadnici a leží na tečně k M0.
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
6
x
y
M0
P
P
M1
x1x0 x0 + x1
2
P má souřadnice:


x0 +
h0
2
, y0 +
h0
2
· f(x0, y0)
y0(x0)




Tečna v bodě P má směrnici:
y (x0 +
h0
2
), tj.
y (x0 +
h0
2
)=
f(x0 +
h0
2
, y0 +
h0
2
·
k1
f(x0, y0)).
Stejnou směrnici by však měla mít i tětiva M0M1 ⇒ souřadnice bodu M1 jsou:
x1 = x0 + h0
y1 = y0 + h0 ·
k2
y (x0 +
h0
2
)
Tyto vztahy lze přepsat do tvaru (obecně)
k1 = f(xn, yn)
k2 = f(xn +
hn
2
, yn +
hn
2
· k1)
yn+1 = yn + hn · k2
Této metodě se říká
modifikovaná Eulerova metoda.
9
Nyní použijeme vlastnost 2)
Známe souřadnice bodu M0. Protože neznáme y-souřadnici bodu M1, nahradíme ho
bodem M1, který má stejnou x-souřadnici a leží na tečně procházející bodem M0.
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
6
7
8
9
x
y
M0
M1
M1
x1x0
M1 má souřadnice:



x0 + h0
=x1
, y0 + h0 ·
ozn.
= k1
f(x0, y0)
=y (x0)




Směrnice tečny v M1 je:
ozn.
= k2
f(x0 + h0, y0 + h0 · f(x0, y0))
Bod M1 dostaneme z podmínky, že směrnice tětivy M0M1 je aritmetickým průměrem
směrnic tečen v M0 a M1, tj.
M1 má souřadnice:
x1 = x0 + h0
y1 = y0 + h0 ·
1
2
(k1 + k2)
Obecně:
k1 = f(xn, yn)
k2 = f(xn + hn, yn + hn · k1)
yn+1 = yn + hn ·
(k1 + k2)
2
Této metodě se říká
Heunova metoda
Poznámka: Obě tyto metody jsou 2.řádu (aproximovali jsme parabolou).
Poznámka: Nejvíce se používá tzv. klasická Runge-Kuttova metoda, která je 4. řádu.
10
Numerické metody pro řešení počátečních
úloh pro ODR – vícekrokové metody
Myšlenka: V jednokrokových metodách se yn+1 počítá pouze s využitím yn
(a hodnot xn, hn). Je rozumné počítat yn+1 s využitím více předchozích hodnot
yn, yn−1, yn−2, . . . , yn−k+1, dosáhneme tím větší přesnosti.
Pro jednoduchost se omezíme na metody s konstantním krokem h (hn = h, ∀ n).
Poznámka: Je třeba si uvědomit, že si lze vymyslet nepřeberné množství metod.
• Jedna z možností je použít metody numerického derivování (špatně podmíněné).
• Další z možností je použít metody numerické integrace
Rovnici y = f(x, y) zintegrujeme od xn do xn+1:
y(xn+1) − y(xn) =
xn+1
xn
f(x, y(x))
=F(x)
dx(5)
Je zřejmé, že funkci F(x) = f(x, y(x)) neznáme. Známe-li ale hodnoty y v bodech
x0, x1, . . . , xn, můžeme vypočítat numerické hodnoty:
F0 = F(x0) = f(x0, y(x0))
F1 = F(x1) = f(x1, y(x1))
...
Fn = F(xn) = f(xn, y(xn))
Pomocí těchto hodnot lze interpolovat funkci F(x) funkcí P(x) a integrál v (5) nahradit
xn+1
xn
P(x) dx.
Interpolace, extrapolace funkce F(x) (2 postupy):
1) F(x) můžeme extrapolovat na intervalu xn, xn+1 pomocí hodnot
F0, F1, . . . , Fn ⇒ explicitně dostaneme y(xn+1) = . . ..
2) F(x) můžeme interpolovat pomocí hodnot F0, F1, . . . , Fn a Fn+1 = F(xn+1) =
f(xn+1, y(xn+1)) ⇒ ve výpočtu integrálu vystoupí yn+1 = y(xn+1) a dostaneme
tak implicitní rovnici s neznámou na obou stranách, tuto rovnici řešíme
postupnými aproximacemi.
11
Adams-Bashfortovy metody
Poznámka: Metody získáné postupem 1).
Postup: Vezmeme posledních k hodnot Fn, Fn−1, . . . , Fn−k+1 a sestrojíme Pk−1(x) interpolační
polynom (k − 1) stupně. Tímto polynomem potom aproximujeme funkci
f(x, y(x)) na intervalu xn, xn+1 , tj. počítáme:
yn+1 = yn +
xn+1
xn
Pk−1(x) dx.
Příklad: Odvoďte vzorec Adams-Bashfortovy metody pro k = 2.
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
6
7
x
F
Fn−1
Fn
xn+1xnxn−1
P1(x)
P1(x) můžeme vyjádřit například
pomocí Langrangeova interpolačního
polynomu:
P1(x) = Fn−1ln−1(x) + Fnln(x),
ln−1(x) =
x − xn
xn−1 − xn
−h
= −
1
h
(x − xn)
ln(x) =
x − xn−1
xn − xn−1
h
=
1
h
(x − xn−1)
P1(x) = Fn−1 −
1
h
(x − xn) + Fn
1
h
(x − xn−1) =
=
1
h
x(Fn − Fn−1) + xnFn−1 − xn−1Fn
xn+1
xn
P1(x) dx =
1
h
(Fn − Fn−1)
x2
n+1
2
−
x2
n
2
1
2 ((xn+h)2
−x2
n)
+(Fn−1xn − Fnxn−1)(xn+1 − xn
h
) =
=
1
h
(Fn − Fn−1)(xnh +
h2
2
) + Fn−1xn − Fnxn−1 =
= Fnxn − Fn−1xn +
h
2
Fn −
h
2
Fn−1 + Fn−1xn − Fnxn−1 =
= Fn(xn − xn−1
h
) +
h
2
Fn −
h
2
Fn−1 = h
3
2
Fn −
1
2
Fn−1 .
⇒ yn+1 = yn +
h
2
(3Fn − Fn−1)
12
Poznámka: Samozřejmě potřebujeme znát prvních k hodnot Fi.
(Ty můžeme vypočítat nějakou jednokrokovou metodou).
Poznámka: Podobně bychom mohli odvodit vzorec Adams-Bashfortovy metody
pro k = 3.
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
x
F
Fn−2 Fn−1
Fn
xn xn+1xn−1xn−2
P2(x)
Opět bychom museli najít interpolační
polynom P2(x) (2. stupně)
a poté zintegrovat přes xn, xn+1 .
Výsledkem je (dcv.):
yn+1 = yn +
h
12
(23Fn − 16Fn−1 + 5Fn−2)
13
Adams-Moultonovy metody
Poznámka: Metody získáné postupem 2).
Postup: Vezmeme posledních k hodnot a přidáme ještě neznámou Fn+1, tj.
Fn+1, Fn, Fn−1, . . . , Fn−k+1. Sestrojíme Qk(x) interpolační polynom k-tého stupně.
Tímto polynomem aproximujeme funkci f(x, y(x)) na intervalu xn, xn+1 , tj. počí-
táme:
yn+1 = yn +
xn+1
xn
Qk(x) dx.
Příklad: Odvoďte vzorec Adams-Moultonovy metody pro k = 1.
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
6
7
x
F
Fn
Fn+1
jako bychom
ji znali
xn+1xn
Q1(x)
Q1(x) můžeme vyjádřit opět např.
pomocí Lagrangeova interpolačního
polynomu:
Q1(x) = Fn+1ln+1(x) + Fnln(x)
ln+1(x) =
x − xn
xn+1 − xn
=
1
h
(x − xn)
ln(x) =
x − xn+1
xn − xn+1
= −
1
h
(x−xn+1)
Q1(x) = Fn+1
1
h
(x − xn) + Fn −
1
h
(x − xn+1) =
=
1
h
[x(Fn+1 − Fn) + Fnxn+1 − Fn+1xn]
xn+1
xn
Qk(x) dx =
1
h
x2
n+1
2
−
x2
n
2
1
2
(xn+1−xn)
h
(xn+1+xn)
(Fn+1 − Fn) + (xn+1 − xn)
h
(Fnxn+1 − Fn+1xn) =
=
1
2
(xn+1 + xn)(Fn+1 − Fn) + Fnxn+1 − Fn+1xn =
=
1
2
xn+1Fn+1 −
1
2
xn+1Fn +
1
2
xnFn+1 −
1
2
xnFn + xn+1Fn − Fn+1xn =
= Fn+1
xn+1
2
+
xn
2
− xn + Fn xn+1 −
xn+1
2
−
xn
2
=
=
h
2
(Fn+1 + Fn)
⇒ yn+1 = yn +
h
2
(Fn+1 + Fn) , kde Fn+1 = f(xn+1, yn+1).
Pozor!
yn+1 = yn +
h
2
f(xn+1, yn+1) + Fn .
Tuto rovnici řešíme iterační metodou např. metodou prosté iterace a tak dostaneme yn+1.
14
Poznámka: Podobně můžeme odvodit vzorec např. pro k = 2
−1 0 1 2 3 4 5 6 7 8
−1
0
1
2
3
4
5
x
F
Fn−1
Fn
Fn+1
xn+1xnxn−1
Q2(x)
Opět bychom museli najít interpolační
polynom Q2(x) (2. stupně).
Poté integrovat přes xn, xn+1
a dostat (dcv.)
yn+1 = yn +
h
12
5 Fn+1 +8Fn − Fn−1
Fn+1=f(xn+1,yn+1)
yn+1 = yn +
h
12
5f(xn+1, yn+1) + 8Fn − Fn−1
Opět vyřešíme iterační metodou → yn+1.
15
Algoritmus prediktor-korektor
Poznámka: Jde o obecné schéma výpočtu.
Princip: Předpokládejme, že máme dostatečně přesně vypočítány hodnoty y0, y1, . . . , yk−1
nějakou explicitní metodou.
Nyní chceme počítat yk.
1) nejprve nějakou explicitní metodou určíme nultou iteraci y
[0]
k jako vstupní hodnotu
pro další výpočet (PREDIKTOR).
2) vypočteme hodnotu pravé strany F
[s]
k = f(xk, y
[s]
k ).
3) vypočteme lepší aproximaci y
[s+1]
k pomocí nějaké implicitní metody s využitím
F
[s]
k =: fk (KOREKTOR).
Pomocí kroků 2) a 3) určíme N iterací y
[1]
k , y
[2]
k , . . . , y
[N]
k (N – dáno).
Na závěr přiřadíme yk = y
[N]
k .
Stejný postup opakujeme pro yk+1, yk+2, . . ..
Poznámka: Dané schéma lze použít na různé metody. Je žádoucí použít explicitní a explicitní
metodu stejného řádu (pro zachování přesnosti). Volba konkrétních metod je
na nás.
Poznámka: Označíme-li operaci:
a) P . . . prediktor
b) E . . . vyčíslení (evaluation)
c) C . . . korektor
Můžeme toto schéma zapsat ve tvaru:
P(EC)N
případně P(EC)N
E, vyčíslujeme-li ještě Fk = f(xk, y
[N]
k ) (což je lepší).
Dostaneme pak různé varianty tohoto schématu:
PEC , PECE
P(EC)2
, P(EC)2
E
P(EC)3
, P(EC)3
E
... ,
...
16
Příklad: Řešte algoritmem prediktor-korektor založeném na Adamsových metodách
druhého řádu na intervalu 0; 0, 6 počáteční úlohu:
y = y + ex
, tj. f(x, y(x)) = y + ex
y(0) = −1
Přesné řešení: y = ex
(x − 1).
Použijeme algoritmus typu PEC.
Vzorec prediktoru má tvar:
y
[0]
n+1 = yn +
h
2
(3Fn − Fn−1)
Korektor:
yn+1 = yn +
h
2
(F
[0]
n+1 + Fn)
Volte krok h = 0, 2.
n xn
přesné
y(xn) y[0]
n F[0]
n yn en
0 0 −1 ••
0 −1← 0
1 0,2 −0,9771 ••
0,2425 •
− 0,9789← 0,0018
2 0,4 −0,8950 P
− 0,9061 E
0,5857→ C
− 0,8960→ 0,0010
3 0, 6 −0,7288 P
− 0,7445 E
1,0776→ C
− 0,7296→ 0,0008
•
Pro určení hodnoty y1 použijeme např. jednokrokovou
modifikovanou Eulerovu metodu (2. řádu):
k1 = f(x0, y0) = y0 + ex0
=
= −1 + 1 = 0
k2 = f(x0 + h/2, y0 + h/2 · k1) =
= −1 + e0,1 .
= 0,1051
y1 = y0 + h · k2
.
=
.
= −1 + 0,2 · 0,1051 = −0,9789
••
Určíme hodnoty F0 a F1.
17