F8370 Moderní metody modelování ve fyzice jaro 2021
D. Hemzal, F. Múnz hemzal@physics.muni.cz
Numerický výpočet derivace
Taylorův rozvoj spojité funkce y{x) v okolí bodu x = a můžeme zapsat jako
y(x) = y(a) + y'(a)(x - a) + -y"(a)(x - a)2 + . . . JJT^^ ~ *f+\ (1)
zde se zbytkem v Lagrangeově tvaru, kde £ leží mezi x a a. Význam zbytku spočívá v tom, že s jeho využitím je rozvoj funkce neaproximativní. Není sice jasné, jak konkrétně zvolit hodnotu £, ale i tak můžeme zbytek prostřednictvím jeho majorizace dobře použít k odhadu přesnosti rozvoje.
Diskretizace
Diskretizované hodnoty nezávislé proměnné x označíme x[n]. V těchto uzlových bodech evaluujeme závislou proměnnou y, označíme přirozeně y(x[n]) = y[n]. V rámci (ekvidistantní) diskretizace zvolme x[n+l]— x[n] = e > 0.
Předpokládejme v dalším speciálně rozvoje jen v uzlových bodech, například
y{x[n] + e) = y{x[n}) + y'{x[n])e + ^y"{x[n])e2 + ..., neboli, s využitím zavedeného značení,
y[n+l] = y[n] + y'[n]e + ^y"[n]e2 + ... (2) Po shrnutí členů s druhou a vyšší derivacemi y do 0(e2) oc e2 můžeme napsat
y'[n] = y[n+1]-y[n]+0(e). (3)
Tomuto vzorci se říká jednobodová dopředná první derivace, podle počtu a typu bodů, které jsou ve vzorci zahrnuty kromě bodu výchozího; chyba této formule je řádu e. Stejně tak jsme ovšem mohli napsat
y{x[n] - e) = y[n] - y'[n]e + ^y"[n]e2 + ... (4)
a byli bychom dostali
y'[n] = y[n]-y[n-1]+0(e), (5)
jednobodovou zpětnou první derivaci. Podotkněme, že numerické hodnoty dopředně a zpětné derivace se obecně liší (a to v řádu chyby 0{e)) i pro spojitou výchozí funkci.
Chybu můžeme zmenšit, když použijeme oba rozvoje současně: přidáním druhé rovnice jsme získali jeden stupeň volnosti navíc a ze soustavy (2), (4) tak můžeme eliminovat člen s druhou derivací,
y[n+l] - y[n-l] 2
y N =-2i-+ ^ ); ( >
tomuto vzorci se říká dvoubodová centrální první derivace. Jak je vidět z tvaru zbytku, pokud bude e malé, skutečně došlo oproti jednobodovým formulím ke zlepšení přesnosti.
Obdobným způsobem se hledají n-bodové k-té derivace. Například dvoubodová centrální druhá derivace se hledá opět ze soustavy (2), (4), tentokrát vyloučíme člen s první derivací a dostáváme
y"[n] = y{^]-2y[n]+y[n-l] + ^ (7)
1
Centrální vzorec nemusí být možné použít u hraničních bodů intervalu, na kterém máme funkci vyjádřenu. Jedna varianta je vrátit se zpět k příslušně orientované derivaci jednobodové, ale (z hlediska chyby) je lepší v takovém případě - například u levé hranice - vypomoci si rozvojem
y[n+2] = y[n] + y'[n]2e + -y"[n]Ac2 + 0{e3)
a ze soustavy (2), (8) opět vyloučit druhou derivaci; obdržíme
y'[n]
3y[n] +4y[n+l] - y[n+2] 2c
Tento vzorec se nazývá dvoubodová dopředná první derivace.
Používání necentrálních vztahů není naštěstí tak časté, jak by se mohlo zdát - zhusta máme na okrajích intervalu definovány okrajové podmínky, které zadávají přímo hodnotu neznámé funkce. Potom sestavujeme derivace již jen ve vnitřních bodech a tam lze užít vzorce centrované. Pokud je předepsána nulová hodnota derivace na hranici, s centrální derivací rovněž vystačíme.
Centrální derivace nižších řádů v ID shrnuje následující obrázek:
Alternativní formule
Pokusme se nyní použít diferenční formule sofistikovanějším způsobem. Předpokládejme, že bychom potřebovali vyjádřit druhé derivace nikoliv pomocí nejbližších sousedů, ale až druhých nejbližších. To je samozřejmě možné, stačí napsat
Výsledek není překvapivý, můžeme si představit, že počítáme diferenci s nejbližšími sousedy ve dvakrát řidší síti. Otázkou nyní zůstává přesnost nového výpočtu, resp. její náprava. Chceme-li dosáhnout i v řidším vzorkování původní chyby, musíme si vypomoci vícebodovou formulí. Napíšeme
y[n±2] = y[n] ± y'[n\2e + -y"[n\Ae2 ± -y"'[n]8e3 + — yIV[n]16£4 + 0((2e)5)
2 o 24
a ze součtu obou vztahů můžeme dostat
y[n+2] -2y[n] +y[n-2]
y[n±A] = y[n] ± y'[n\As + -y"[n]16e2 ± -y"'[n]64e3 + —yIV[n]256£4 + 0((Ae)5)
2 b 24
odkud (s majorizací chyby) můžeme napsat
y"[n]
y[n+A] + 16y[n+2] - 30y[n] + 16y[n-2] - y[n-A] 48^2
+ 0((Aef)
(tedy opět běžná čtyřbodová formule v řidší síti).
Více dimenzí
Ve více dimenzích musíme použít odpovídající tvar Taylorova rozvoje:
/(x) = /(xq) + d/(x0) + -d2/(x0) + • • •
kde
d f = ^ • Ax = fxAx + fyAy + ...
2
a diferenciály vyšších řádů se získají formálním umocněním. Ve dvou dimenzích s (ekvidistantními) kroky sítě ex, £y tedy například platí
/(x) = /(xo) + fx{*o)£x + fy{*o)£y + \ {fxxi^O^l + 2fxy(x0)cxCy + fyy(xo)£y) + • • •
Tuto formuli můžeme využít k vyjádření dvourozměrného Laplasiánu jako zobecnění druhé derivace. Zaměříme-li se na pravoúhlý kříž kolem uzlu a;[/c, Z], můžeme ve zjednodušeném případě homogenní kartézské sítě {ex = By = e) psát
f[k+l,l] = f[k, Z] + fx[k, l]e + \fxx[k, l]s2 + 0(£3) f[k-l,l] = f[k, Z] - fx[k, l]e + \fxx[k, l]s2 + 0(£3) f[k,l+l] = f[k, 1} + fy[k, l]e + l-fyy[k, l]s2 + 0(£3)
f[k,l~l]   =   f[k, l]   -   fy[k, l]8 +   -fyy[k, l]^   + O^)
Sečtením všech rovnic dostáváme hledanou aproximaci
A/[M] =fxx[k,l]+fyy[k,l]
f[k + l,l] + f[k, l + 1] - 4f[k,l] + f[k-l,l] + f[k,l-l]
+ 0(e)
Ze struktury výpočtu je patrné, že jsme nejprve mohli sečíst první dvě rovnice a získat tak známou formuli pro druhou derivaci v ID, následně provést totéž ve druhé dimenzi a oba mezivýsledky složit.
Uvedená aproximace Laplasiánu je kanonická, avšak zdaleka ne jediná. Předpokládejme, že bychom se místo pravoúhlého kříže rozhodli použít kříž ůhlopříčný. Potom bychom mohli psát
/[fc+l,Z+l] = f[k, 1} + fx[k, l]e + fy[k, l]e + \ (fxx[k, 1} + 2fxy[k, 1} + fyy[k, l])s2 + 0(s3)
f[k-l,l+l] = f[k, 1} - fx[k, l]e + fy[k, l]e + - (fxx[k, 1} - 2fxy[k, 1} + fyy[k, l])e2 + 0(s3)
f[k-l,l-l] = f[k, 1} - fx[k, l]e - fy[k, l]e + \ (fxx[k, 1} + 2fxy[k, 1} + fyy[k, l]) s2 + 0(£3)
f[k+l,l-l] = f[k, 1} + fx[k, l]e - fy[k, l]e + \ (fxx[k, 1} - 2fxy[k, 1} + fyy[k, l])s2 + 0(£3)
Nyní opět pouhým sečtením všech rovnic dostaneme
/[fc+l,ž+l] +/[fc-l,ž+l] -4/[fc,ž] +/[fc-l,ž-l] +/[fc+l,ž-l]        f e A/[fc, l] =-—2- V2
a tato formule je (mimo jiné) o něco přesnější než předchozí.
Grafické znázornění obou získaných aproximací 2D Laplasiánu spolu s příslušnou smíšenou derivací zachycuje následující obrázek:
82A
2e2A
4e-
2 ď2 dxdy
3
Pokud bychom pokračovali s výpočtem kanonického Laplasiánu do více dimenzí, obdržíme posloupnost
Adaptace sítě na symetrii úlohy
Konečné diference lze poměrně dobře využít i ve složitějších souřadných systémech než kartézských. Následující obrázky zachycují výpočty Laplasiánu ve voštinové síti a v polárních souřadnicích:
4