Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Drsná matematika III – 5. přednáška
Diferenciální operátory
Jan Slovák
Masarykova univerzita
Fakulta informatiky
17. 10. 2011
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Obsah přednášky
1 Literatura
2 Modely založené na změnách
Lineární a nelineární modely
3 Obyčejné diferenciální rovnice
Rovnice se separovanými proměnnými
Systémy ODE prvního řádu
Lineární diferenciální rovnice
Tlumený oscilátor
Lineární diferenciální rovnice s konstantními koeficienty
4 Numerické metody
Eulerova metoda
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Kde je dobré číst?
Riley, K.F., Hobson, M.P., Bence, S.J. Mathematical Methods
for Physics and Engineering, second edition, Cambridge
University Press, Cambridge 2004, ISBN 0 521 89067 5, xxiii
+ 1232 pp.
Kurz 18.03 na MIT OpenCourseWare, viz
http//ocw.mit.edu/OcwWeb/Mathematics
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
lineární model prvního řádu
Derivace pracuje s okamžitými změnami studovaných veličin. Ze
stejných důvodů jsme kdysi zaváděli diference a právě vztahy mezi
hodnotami veličin a změnami těch samých nebo jiných veličin vedly
k rovnicím.
Nejjednodušším modelem bylo úročení vkladů nebo půjček a totéž
pro tzv. Malthusiánský model populace. Přírůstek byl úměrný
hodnotě.
V rámci spojitého modelování by stejný požadavek vedl na rovnici
vztahující derivaci funkce y (x) s její hodnotou
y (x) = r · y(x)
s konstantou úměrnosti r. Je snadné uhodnout řešení této rovnosti
y(x) = C erx
s libovolnou konstantou C.
Tuto konstantu určíme jednoznačně volbou tzv. počáteční
hodnoty y0 = y(x0) v nějakém bodě x0.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Pokud je část růstu v našem modelu dána konstatním působením
nezávislém na hodnotě y nebo x (např. paušální poplatky za vedení
účtu nebo přirozený úbytek populace třeba v důsledku porážek na
jatkách), přidáme na pravé straně konstantu s:
y (x) = r · y(x) + s.
Zjevně bude řešením této rovnice funkce
y(x) = C erx
−
s
r
.
K tomuto závěru je velice lehké dojít, pokud si uvědomíme, že
množinou všech řešení první (homogenní) rovnice je jednorozměrný
vektorový prostor, zatímco řešení druhé (nehomogenní) rovnice se
obdrží přičtením kteréhokoliv jednoho jejího řešení ke všem řešením
předchozí rovnice. Lze pak snadno najít konstantní řešení y(x) = k
pro k = −s
r .
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Nelineární model prvního řádu
U diferenčních rovnic jsme diskutovali tzv. logistický model
populačního růstu založený na předpokladu, že poměr změny
velikosti populace p(n + 1) − p(n) a její velikosti p(n) je v afinní
závislosti na samotné velikosti populace.
Nyní bychom tentýž vztah pro spojitý model patrně formulovali pro
populaci p(t) závislou na čase t jako
p (t) = p(t) −
r
K
p(t) + r ,
tj. při hodnotě p(t) = K pro velkou konstantu K je přírůstek
nulový, zatímco pro p(t) blízké nule je poměr rychlosti růstu
populace k její velikosti blízký r (což je malé číslo v řádu setin
vyjadřující rychlost růstu populace za dobrých podmínek.)
Derivováním lze přímo ověřit, že následující funkce je řešením pro
každou konstantu C
p(t) =
K
1 + CK e−rt
.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Srovnání diskrétního a spojitého modelu
y
100
80
60
40
20
t
0
200150100500
100
80
60
40
20
x
200150100500
Srovnáním červeného grafu řešení s K = 100, r = 0, 05 a C = 1
(volba C odpovídá p(0) = 1) s řešením diferenční rovnice (napravo)
vidíme, že skutečně oba přístupy k modelování populací dávají
docela podobné výsledky. Pro srovnání je do levého obrázku zeleně
vkresleno řešení Malthusiánského modelu s odpovídajícími daty.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Definition
Diferenciální rovnice prvního řádu je
F(y (x), y(x), x) = 0
s nějakou pevnou funkcí F, která každé trojici reálných čísel přiřadí
jedno reálné číslo.
Zápis připomíná implicitně zadané funkce y(x), nicméně navíc je tu
závislost na derivaci hledané funkce y(x).
Obvykle pracujeme s rovnicemi, které jsou vyřešeny vzhledem k
derivaci, tj.
y (x) = f (x, y(x)),
dále se omezíme jen na tento případ.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Taková rovnice zadává pro každou hodnotu (x, y) v rovině vektor
(1, f (x, y)), tj. rychlost se kterou nám rovnice grafu řešení přikazuje
pohybovat se rovinou. Např. pro logistický model
p (t) = p(t) −
r
K
p(t) + r
0
x
200150100500
y(x)
100
80
60
40
20
Intuitivně lze na základě takových obrázků očekávat, že pro každou
počáteční podmínku bude existovat právě jedno řešení naší rovnice.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Existence i jednoznačnost řešení skutečně platí pro všechny
rozumné funkce f , my si výsledek sformulujeme pro dosti velkou
třídu rovnic takto:
Theorem (O existenci a jednoznačnosti řešení ODE)
Nechť funkce f (x, y) : R2 → R má spojité parciální derivace. Pak
pro každý bod (x0, y0) ∈ R2 existuje interval [x0 − a, x0 + a], s
a ∈ R kladným, a právě jedna funkce y(x) : R → R, která je
řešením rovnice
y (x) = f (x, y(x)).
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Náznak důkazu
Funkce y(x) je řešením naší rovnice tehdy a jen tehdy, když
y(x) = y0 +
x
x0
y (x) dx = y0 +
x
x0
f (x, y(x)) dx.
Pravá strana tohoto výrazu je ovšem, až na konstantu, integrální
operátor
L(y)(x) = y0 +
x
x0
f (x, y(x)) dx
a při řešení diferenciální rovnice hledáme pevný bod pro tento
operátor L, tj. chceme najít funkci y = y(x) s L(y) = y.
Důkaz spočívá v odhadu, že pro dostatečně malý interval kolem x0
bude takový operátor zmenšovat vzdálenost funkcí. Z obecné věty o
kontrakci pak vyplývá hledané tvrzení.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Užitečným typem rovnic, pro který známe elementární postup k
řešení jsou tzv. rovnice se separovanými proměnnými:
y (x) = f (x) · g(y(x))
pro dvě dostatečně hladké funkce jedné reálné proměnné f a g.
Obecné řešení tu lze získat integrací, tj. nalezením primitivních
funkcí
G(y) =
dy
g(y)
, F(x) = f (x)dx.
Pak totiž spočtením funkce y(x) z implicitně zadaného vztahu
F(x) + C = G(y) s libovolnou konstantou C vede k řešení.
Skutečně, derivováním této rovnosti (s použitím pravidla pro
derivování složené funkce G(y(x)) dostaneme skutečně
1
g(y) · y (x) = f (x).
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Příklad rovnice
Najdeme řešení rovnice
y (x) = x · y(x).
Přímým výpočtem dostaneme ln |y(x)| = 1
2x2 + C. Odtud to
vypadá (alespoň pro kladná y) na
y(x) = e
1
2
x2+C
= D · e
1
2
x2
,
kde D je nyní libovolná kladná konstanta.
Ve skutečnosti dostaneme všechna řešení, když uvážíme všechna
reálná D
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
y(x)
x
3
3
2
1
2
0
-1
1
-2
-3
0-1-2-3
y(x)
x
3
3
2
1
2
0
-1
1
-2
-3
0-1-2-3
Na obrázku jsou vynesena dvě řešení, která ukazují na nestabilitu
rovnice vůči počátečním podmínkám: Jestliže pro libovolné x0
volíme y0 blízké nule, pak se nám dramaticky mění chování
výsledného řešení. Navíc si povšimněme konstatního řešení
y(x) = 0, které odpovídá počáteční podmínce y(x0) = 0.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Jestliže lehce pozměníme rovnici na
y (x) = 1 − x · y(x),
narazíme naopak na stabilní chování viditelné na následujícím
obrázku. Tuto rovnici už ale neumíme řešit pomocí separace
proměnných.
y(x)
x
3
6
2,5
2
5
1,5
1
4
0,5
0
3210
y(x)
x
3
6
2,5
2
5
1,5
1
4
0,5
0
3210
Zato umíme stejným postupem řešit logistický model populace.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Na řešení rovnice y (x) = f (x, y) lze také pohlížet jako na hledání
(parametrizované) křivky (x(t), y(t)) v rovině, kde jsme již předem
pevně zvolili parametrizaci proměnné x(t) = t. Takto můžeme
nejen zapomenout na tuto pevnou volbu pro jednu proměnnou x,
nýbrž hlavně přibrat libovolný počet proměnných.
Například v rovině můžeme psát takový systém ve tvaru
x (t) = f (t, x(t), y(t)), y (t) = g(t, x(t), y(t))
se dvěmi funkcemi f , g : R3 → R se spojitými derivacemi. Obdobně
pro více proměnných.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Jednoduchým příkladem v rovině může sloužit systém rovnic
x (t) = −y(t), y (t) = x(t).
Snadno lze uhádnout (nebo aspoň ověřit), že řešením takového
systému je
x(t) = R cos t, y(t) = R sin t
s libovolnou nezápornou konstantou R a křivky řešení budou právě
parametrizované kružnice o poloměru R.
Na takové systémy umíme přímo rozšířit platnost věty o
jednoznačnosti a řešení:
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Theorem (O existenci a jednoznačnosti řešení systémů ODE)
Nechť funkce fi (t, x1, . . . , xn) : Rn+1 → R, i = 1, . . . , n všechny
mají spojité parciální derivace. Pak pro každý bod
(t0, z1, . . . , zn) ∈ R2 existuje interval [t0 − a, t0 + a], s a ∈ R
kladným, a právě jedna funkce y(t) : R → Rn, která je řešením
systému rovnic
x1(x) = f1(t, x1(t), . . . , xn(x)), . . . , xn(x) = fn(t, x1(t), . . . , xn(x))
s počáteční podmínkou
x1(t0) = z1, . . . , xn(t0) = zn.
Ve skutečnosti je možné se omezit na tzv. autonomní systémy
rovnic, tj. ty kde pravá strana není explicitně závislá na čase t.
Obecný případ ve větě se z nich dostane stejným způsobem, jako
jsme pracovali v rovině s jednou rovnicí y (x) = f (x, y).
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Poslední systém dvou autnomních rovnic, jehož řešení jsou
parametrizované kružnice je dobrým příkladem.
Samotný systém rovnic (a každý obecný také) si pak můžeme
představit jako „pole vektorů“ v rovině, prostoru atd., které udává
„tok prostředí“. Řešením je pak křivka, která odpovídá skutečnému
toku jedné „částečky“ tohoto prostředí.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Model „Lotka – Voltera“ pro dravce a kořist
Jako o něco složitější příklad si uveďme klasický populační model
„dravec – kořist“, který zavedli ve dvacátých létech minulého století
pánové Lotka a Volterra.
Označíme x(t) vývoj počtu jedinců v populaci kořisti a y(t) totéž
pro dravce. Přepokládáme, že přírůstek kořisti by se řídil
Malthusiánským modelem (tj. exponenciální růst), kdyby nebyli
loveni. U dravce naopak očekáváme, že by bez kořisti pouze
přirozeně vymíral (tj. exponenciální pokles stavů). Přitom ale ještě
musíme uvážit interakci dravce s kořistí, kterou očekáváme přímo
úměrnou počtu obou. Dostáváme tak tzv. Lotka–Volterra model
x (t) = αx(t) − βy(t)x(t)
y (t) = −γy(t) + δx(t)y(t)
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Operace derivování je lineární zobrazení z (dostatečně) hladkých
funkcí do funkcí. Pokud derivace ( d
dx )j jednotlivých řádů j
vynásobíme pevnými funkcemi aj (x) a výrazy sečteme, dostaneme
tzv. lineární diferenciální operátor:
y(x) → D(y)(x) = ak(x)y(k)
(x) + · · · + a1(x)y (x) + a0y(x).
Řešit příslušnou homogenní lineární diferenciální rovnici pak
znamená najít funkci y splňující D(y) = 0, tj. obrazem je identicky
nulová funkce.
Ze samotné definice je zřejmé, že součet dvou řešení bude opět
řešením, protože pro libovolné funkce y1 a y2 platí
D(y1 + y2)(x) = D(y1)(x) + D(y2)(x).
Obdobně je také konstantní násobek řešení opět řešením. Celá
mmnožina všech řešení lineární diferenciální rovnice k-tého řádu je
tedy vektorovým prostorem.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Přímou aplikací předchozí věty o jednoznačnosti a existenci řešení
rovnic dostáváme:
Corollary
Vektorový prostor všech řešení homogenní lineární diferenciální
rovnice k–tého řádu je vždy dimenze k. Proto můžeme vždy řešení
zadat jako lineární kombinaci libovolné množiny k lineárně
nezávislých řešení. Taková řešení jsou zadána jednoznačně lineárně
nezávislými počátečními podmínkami na hodnotu funkce y(x) jejích
prvních (k − 1) derivací.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Zkusme si popsat jednoduchý model pro pohyb nějakého tělesa
upnutého k jednomu bodu silnou pružinou. Je-li y(t) výchylka
našeho tělesa od bodu y0 = y(0) = 0, pak lze uvažovat, že
zrychlení y (t) v čase t bude úměrné velikosti výchylky, avšak s
opačným znaménkem. Dostáváme tedy tzv. rovnici oscilátoru
y (t) = −y(t).
Tato rovnice odpovídá systému rovnic
x (t) = −y(t), y (t) = x(t).
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Řešením takového systému je
x(t) = R cos(t − τ), y(t) = R sin(t − τ)
s libovolnou nezápornou konstantou R, která určuje maximální
amplitudu, a konstantou τ, která určuje fázový posun.
Pro určení jednoznačného řešení potřebujeme proto znát nejen
počáteční polohu y0, nýbrž také rychlost pohybu v tomto okamžiku.
Těmito dvěma údaji bude určena jak amplituda tak fázový posun
jednoznačně.
Předpokládejme, že vlivem vlastností materiálu pružiny bude také
působit síla úměrná okamžité rychlosti pohybu našeho objektu, s
opačným znaménkem neš amplituda.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
To vyjádříme dodatečným členem s první derivací a naše rovnice je
y (t) = −y(t) − αy (t),
kde α je konstanta, která vyjadřuje velikost tlumení. Na obrázku
jsou fázové diagramy pro řešení s dvěmi různými počátečními
podmínkami. Nalevo je nulové tlumení, napravo je α = 0.3
0
5
-3-3
10-2
-2
t-1
-1
0 15
0
1x(t)
y(t)
1
2
203
2
3
Tlumeneoscilace
0
5
-3-3
10-2
-2
t-1
-1
0 15
0
1x(t)
y(t)
1
2
203
2
3
Tlumeneoscilace
Samotné oscilace jsou vyjádřeny hodnotami na ose y, hodnoty x
zobrazují rychlost pohybu.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Připomeňme homogenními lineární diferenční rovnice.
Analogie jde i dále v okamžiku, kdy jsou všechny koeficienty aj
diferenciálního operátoru D konstantní. Už jsme viděli u takové
rovnice prvního řádu, že řešením je exponenciála s vhodnou
konstantou u argumentu. Stejně jako u diferenčních rovnic se
podbízí vyzkoušet, zda takový tvar řešení y(x) = eλx s neznámým
parametrem λ může splnit rovnici k–tého řádu. Dosazením
dostaneme
D(eλx
) = akλk
+ ak−1λk−1
+ · · · + a1λ + a0(x) eλx
.
Parametr λ tedy vede na řešení lineární diferenciální rovnice s
konstantními koeficienty tehdy a jen tehdy, když je λ kořenem tzv.
charakteristického polynomu akλk + · · · + a1λ + a0.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Pokud má charakteristický polynom k různých kořenů, dostáváme
bázi celého vektorového prostoru řešení. Pokud je λ násobný kořen,
přímým výpočtem s využitím toho, že je pak také kořenem derivace
charakteristického polynomu, dostaneme, že je řešením i funkce
x eλx . Podobně pak pro vyšší násobnost dostáváme různých
řešení eλ, x eλx , . . . , x eλx .
U obecné lineární diferenciální rovnice předepisujeme nenulovou
hodnotu diferenciálního operátoru D. Opět úplně analogicky k
úvahám o systémech lineárních rovnic nebo u lineárních
diferenčních rovnic přímo vidíme, že obecné řešení takovéto
(nehomogenní) rovnice
D(y)(x) = b(x)
pro nějakou pevně zadanou funkci b(x) je součtem jednoho
jakéhokoliv řešení této rovnice a množiny všech možných řešení
příslušné homogenní rovnice D(y)(x) = 0. Celý prostor řešení je
tedy opět pěkný konečněrozměrný afinní prostor, byť ukrytý v
obrovském prostoru funkcí.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Kromě tak jednoduchých rovnic, jako jsou ty lineární s konstantními
koeficienty se v praxi většinou setkáváme s postupy, jak přibližně
spočíst řešení rovnice, se kterou pracujeme.
Už jsme podobné úvahy dělali všude tam, kde jsme se zabývali
aproximacemi (tj. zejména lze doporučit porovnání s dřivějšími
odstavci o splajnech, Taylorových polynomech a Fourierových
řadách). S trochou odvahy můžeme také považovat diferenční a
diferenciální rovnice za vzájemné aproximace. V jednom směru
nahrazujeme diference diferenciály (např. u ekonomických nebo
populačních modelů), ve druhém pak naopak.
Zastavíme se na chvilku u nahrazování derivací diferencemi.
Nejdříve si však připomeneme obvyklé značení pro zápis odhadů
chyb.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Definition
Pro funkci f (x) v proměnné x řekneme, že je v okolí hromadného
bodu x0 svého definičního oboru řádu velikosti O(ϕ(x)) pro
nějakou funkci ϕ(x), jestliže existuje okolí U bodu x0 a konstanta
C taková, že
|f (x)| ≤ C · |ϕ(x)|
pro všechny x ∈ U. Limitní bod x0 bývá často i nevlastní hodnota
±∞.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Nejobvyklejší příklady jsou O(xp) pro polynomiální řád velikosti
a to v nule nebo v nekonečnu, O(ln x) pro logaritmický řád
velikosti v nekonečnu atd. Všimněme si, že logaritmický řád
velikosti nezávisí na volbě základu.
Dobrým příkladem je aproximace funkce jejím Taylorovým
polynomem řádu k v bodě x0. Taylorova věta pro funkce jedné
proměnné říká, že chyba této aproximace je O(hk+1), kde h je
přírůstek argumentu x − x0 = h.
Podobné úvahy jsme dělali i u Fourierových řad.
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
V případě obyčejných diferenciálních rovnic je nejjednodušším
schématem aproximace tzv. Eulerovými polygony. Budeme ji
prezentovat pro jednu obyčejnou rovnici s jednou nezávislou a
jednou závislou veličinou. Úplně stejně ale funguje pro systémy
rovnic, když skalární veličiny a jejich derivace v čase t nahradíme
vektory závislé na času a jejich derivacemi.
Uvažujme tedy opět rovnici (pro jednoduchost a bez újmy na
obecnosti prvního řádu)
y (t) = f (t, y(t)).
Literatura Modely založené na změnách Obyčejné diferenciální rovnice Numerické metody
Označme si diskrétní přírůstek času h, tj. tn = t0 + nh, a
yn = y(tn). Z Taylorovy věty (se zbytkem druhého řádu) a naší
rovnice vyplývá, že
yn+1 = yn + y (tn)h + O(h2
) = yn + f (tn, yn)h + O(h2
).
Jestliže tedy od t0 do tn uděláme n takových kroků o přírůstek h,
bude očekávaný odhad celkové chyby vyplývající z lokálních
nepřesností naší lineární aproximace nejvýše hO(h2), tj. chyba bude
v řádu velikosti O(h). Ve skutečnosti vstupují při výpočtu do hry
ještě zaokrouhlovací chyby.
Při numerickém řešení Eulerovou metodou postupujeme tak, že za
přibližné řešení považujeme po částech lineární polygon definovaný
výše.