Programování pro fyziky 2021/22 – Numerické recepty č. 5

Náš průvodní text je především v poznámkách k přednášce. Pomáhat si budeme i odkazy do knih Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR), které jsou (v Second Edition) na webu volně k dispozici. Vidět je ovšem lze jen z Adobe Readeru po instalaci pluginu FileOpen. Interní PDF prohlížeče Firefoxu a Chromu nám nepomohou.

Soustavy obyčejných diferenciálních rovnic

(NR kap. 16) Úlohou je nalézt numerickou aproximaci funkcí, pro které je předepsána soustava \(R\) obyčejných diferenciálních rovnic prvního řádu,

\[\frac{dy_r(x)}{dx}=f_r(x,y_1(x),\dots,y_R(x)),\ r=1,\dots,R, \quad\mbox{neboli}\quad \frac{d\mathbf{y}(x)}{dx}=\mathbf{f}(x,\mathbf{y}(x)),\]

s počáteční podmínkou

\[y_r(x_0)=y_{0r} \quad\mbox{neboli}\quad \mathbf{y}(x_0)=\mathbf{y}_0.\]

Soustavy vyššího řádu lze převádět na soustavy prvního řádu, obdobně jako v následující ukázce. Numerickým řešením úlohy na síti \(x_n,\ n=0,1,\dots,\) je pro každé \(r=1,\dots,R\) posloupnost \(y_{nr}\), která lépe či hůře aproximuje přesné řešení \(y_r(x_n)\).

Ukázka: rovnice pro trajektorii \(X(t)=(x(t),z(t))\) bodu o hmotnosti \(M\) v silovém poli \(F=(F_x,F_z)\) se zadanou polohou a rychlostí \(V(t)=(v_x,v_z)\) v počátečním čase \(t_0\),

\[M\frac{d^2X(t)}{dt^2}=F(t,X(t)),\ X(t_0)=X_0,\ V(t_0)=V_0,\]

lze vyjádřit jako soustavu prvního řádu,

\[\begin{split}\frac{d}{dt}\left(\begin{array}{r} x(t) \\ z(t) \\ v_x(t)\\ v_z(t) \\ \end{array}\right) = \left(\begin{array}{l} v_x(t) \\ v_z(t) \\ F_x(t,x,z,v_x,v_z)/M \\ F_z(t,x,z,v_x,v_z)/M \\ \end{array}\right).\end{split}\]

Eulerova metoda

Dvě podoby nejjednodušší numerické metody lze získat použitím nultého a prvního členu Taylorova rozvoje pro \(y(x)\) v okolí \(x_n\), resp. \(x_{n+1}\) (vynecháváme index rovnice \(r\)):

\[y(x_n+h)\approx y(x_n)+hy'(x_n), \quad\mbox{resp.}\quad y(x_{n+1}-h)\approx y(x_{n+1})-hy'(x_{n+1}).\]

Když zde položíme \(x_n+h=x_{n+1}\), nahradíme výrazy \(y(x_n)\) pro přesné řešení odpovídajícími výrazy \(y_n\) pro numerické řešení a dosadíme za \(y'(x)\) pravou stranu výchozích rovnic \(f(x,y_1,\dots,y_R)\), získáme následující:

explicitní Eulerova metoda pro 1 rovnici

\[y_{n+1}=y_n+hf(x_n,y_n),\quad n=0,1,\dots,\]

implicitní Eulerova metoda pro 1 rovnici

\[y_{n+1}=y_n+hf(x_{n+1},y_{n+1}).\]

Pro soustavu \(R\) rovnic mají explicitní a implicitní Eulerova metoda tvar

\[y_{n+1,r}=y_{nr}+hf(x_n,y_{n1},\dots,y_{nR}),\ r=1,\dots,R, \quad\mbox{neboli}\quad \mathbf{y}_{n+1}=\mathbf{y}_n+h\mathbf{f}(x_n,\mathbf{y}_n),\]
\[y_{n+1,r}=y_{nr}+hf(x_{n+1},y_{n+1,1},\dots,y_{n+1,R}), \quad\mbox{neboli}\quad \mathbf{y}_{n+1}=\mathbf{y}_n+h\mathbf{f}(x_{n+1},\mathbf{y}_{n+1}).\]

Přesnost mají v taylorovském smyslu obě Eulerovy metody stejnou (prvního řádu), podstatně lepší je však stabilita implicitní varianty (což platí pro příbuzné explicitní a implicitní metody vcelku obecně). To lze ukázat na úloze

\[y'=-cy,\ y(0)=1,\ c>0,\]

s analytickým řešením \(y(x)=e^{-cx}\) (viz NR kap. 16.6):

  1. explicitní metoda: \(y_{n+1}=y_n+h(-cy_n)=(1-ch)y_n\), což dává nestabilní (divergující) numerické řešení při \(|1-ch|>1\), čili při kroku \(h>2/c\);

  2. implicitní metoda: \(y_{n+1}=y_n+h(-cy_{n+1})=y_n/(1+ch)\), což nediverguje ani pro \(h\to\infty\).

Octave: vývoj \(y_n\) pro \(y'=-cy\) pomocí explicitní (pole ye, red) a implicitní (yi, blue) Eulerovy metody při \(h>2/c\):

c=1; x0=0; y0=1; ye=y0; yi=y0; xmax=10.4; nmax=5; h=xmax/nmax; x=x0:h:xmax;
for n=1:nmax, ye(n+1)=(1-c*h)*ye(n); yi(n+1)=yi(n)/(1+c*h); end; plot(x,exp(-c*x),'g',x,ye,'r',x,yi,'b')

Cenou za lepší stabilitu implicitních metod je jejich větší numerická náročnost. Ta je způsobena buď několika iteracemi v každém kroku, kdy se po inicializaci \(y_{n+1}^{[0]}\) explicitní metodou opakovaně volá implicitní metoda, \(y_{n+1}^{[j+1]}=y_n+hf(x_{n+1},y_{n+1}^{[j]}),\) \(j=0,1,\dots,\) nebo řešením soustavy algebraických rovnic v každém kroku, k čemuž vede linearizace pravé strany soustavy (aplikace Newtonovy metody), zde \(\mathbf{f}=\mathbf{f}(\mathbf{y})\), ve vzorci implicitní metody:

\[\mathbf{y}_{n+1}=\mathbf{y}_n+h\mathbf{f}(\mathbf{y}_{n+1}) =\mathbf{y}_n+h\left(\mathbf{f}(\mathbf{y}_n)+\mathbf{J}\cdot(\mathbf{y}_{n+1}-\mathbf{y}_n)\right),\]

kde Jacobiovu matici \(\frac{\partial\mathbf{f}}{\partial\mathbf{y}}|_{\mathbf{y}_n}\) značíme \(\mathbf{J},\ J_{ij}=\frac{\partial f_i}{\partial y_j}\). Pak

\[(\mathbf{I}-h\mathbf{J})\cdot\mathbf{y}_{n+1}=(\mathbf{I}-h\mathbf{J})\cdot\mathbf{y}_n+h\mathbf{f}(\mathbf{y}_n)\]

a

\[\mathbf{y}_{n+1}=\mathbf{y}_n+h(\mathbf{I}-h\mathbf{J})^{-1}\cdot\mathbf{f}(\mathbf{y}_n).\]

Právě odvozený vzorec se nazývá semi-implicitní Eulerovou metodou.

Rungovy-Kuttovy metody

(NR kap 16.1) Obecný vzorec (zde: explicitních) Rungových-Kuttových (RK) metod (zde: pro 1 rovnici) má tvar

\[\begin{split}\begin{array}{rcl} y_{n+1} &=& y_n+h\sum_{i=1}^p c_ik_i, \\ k_1 &=& f(x_n,y_n), \\ k_i &=& f(x_n+h\alpha_i,y_n+h\sum_{j=1}^{i-1}\beta_{ij}k_j),\ i=2,\dots,p. \end{array}\end{split}\]

Parametry \(c_i,\ \alpha_i,\ \beta_{ij}\) konkrétní metody se volí tak, aby vážený průměr několika směrnic \(k_i\) podél hledané funkce \(y(x)\) na intervalu \(\langle x_n,x_{n+1}\rangle\) aproximoval Taylorův rozvoj \(y(x)\) co nejvyššího řádu. Konstanty \(c_i\) jsou vahami směrnic, tedy \(\sum_{i=1}^p c_i=1\), konstanty \(\alpha_i\) jsou z intervalu \(\langle0,1\rangle\) a konstanty \(\beta_{ij}\) vyplňují dolní trojúhelníkovou matici. Počet \(p\) vyčíslení pravé strany odpovídá pro \(p\le4\) řádu přesnosti metody.

RK1 metoda prvního řádu přesnosti má parametry \(p=1,\ c_1=1\) a je to explicitní Eulerova metoda. Pro každé \(p>1\) existuje \(\infty\) RK metod, používá se jich však jen několik s jednoduše vyhlížejícími parametry nebo speciálními vlastnostmi. Dvě populární metody druhého řádu (NR 16.1.2):

RK2 metoda středního bodu (midpoint method)

\[y_{n+1}=y_n+hf(x_n+\frac12h,y_n+\frac12hf(x_n,y_n))\,,\]

RK2 Heunova metoda (Heun-Verfahren)

\[y_{n+1}=y_n+\frac12h(f(x_n,y_n)+f(x_n+h,y_n+hf(x_n,y_n))\,.\]

Odvození RK2 metod: Taylorův rozvoj \(y(x)\) druhého řádu je \(y_{n+1}\approx y_n+hy'+\frac12h^2y''\), což při \(y'=f(x_n,y_n)\equiv f\) a \(y''=f'=\partial_x f+f\partial_y f\) lze přepsat na tvar \(y_{n+1}\approx y_n+hf+\frac12h^2\partial_xf+\frac12h^2f\partial_yf\). Obecný vzorec RK2 metody \((p=2)\) je \(y_{n+1}=y_n+hc_1f(x_n,y_n)+hc_2f(x_n+h\alpha_2,y_n+h\beta_{21}f)\approx y_n+hc_1f+hc_2(f+h\alpha_2\partial_xf+h\beta_{21}f\partial_yf)\), kde jsme použili Taylorův rozvoj funkce 2 proměnných. Srovnáním členů v obou vzorcích získáme soustavu 3 nelineárních rovnic pro 4 parametry metody: \(c_1+c_2=1\), \(c_2\alpha_2=\frac12\), \(c_2\beta_{21}=\frac12\). Můžeme volit \(c_2=C\), pak \(c_1=1-C\), \(\alpha_2=\beta_{21}=1/(2C)\). Pro \(C=1\) máme metodu středního bodu, pro \(C=\frac12\) Heunovu metodu.

Populární RK4 metoda čtvrtého řádu (NR 16.1.3):

\[p=4,\ c=(\frac16,\frac13,\frac13,\frac16),\ \alpha=(-,\frac12,\frac12,1),\ \beta=(-,-,-,-|\ \frac12,-,-,-|\ 0,\frac12,-,-|\ 0,0,1,-)\]

neboli \(y_{n+1}=y_n+h\left(\frac16k_1+\frac13k_2+\frac13k_3+\frac16k_4\right),\)

kde \(k_1=f(x_n,y_n),\ k_2=f(x_n+\frac12h,y_n+\frac12hk_1),\ k_3=f(x_n+\frac12h,y_n+\frac12hk_2),\ k_4=f(x_n+h,y_n+hk_3).\)

Odvození se vede podobně jako pro RK2, získá se soustava 11 nelineárních rovnic pro 13 parametrů metody.

Pro pravou stranu \(f(x)\) nezávislou na \(y(x)\) je krok Eulerovy metody ekvivalentní (levostrannému) obdélníkovému kvadraturnímu vzorci, \(y_{n+1}=y_n+hf_n\), RK2 Heunova metoda odpovídá lichoběžníkovému pravidlu, \(y_{n+1}=y_n+\frac12h(f_n+f_{n+1})\), a RK4 metoda Simpsonovu pravidlu s krokem \(\frac h2\): \(y_{n+1}=y_n+\frac h6(k_1+2k_2+2k_3+k_4)\), kde \(k_1=f(x_n),\ k_2=k_3=f(x_n+\frac h2),\ k_4=f(x_{n+1})\).

MATLAB: např. funkce ode45(f,tspan,y0) pro RK metodu zapouzdřující vzorce 4. a 5. řádu přesnosti (při \(p=6\)). To umožňuje provádět v každém kroku odhad chyby a adaptivně tak regulovat délku (následujícího nebo revidovaného) kroku podle předepsané mezní chyby (NR kap. 16.2).

Octave: funkce lsode(fcn,y0,xvec) z volně dostupné knihovny ODEPACK, kde xvec je vektor popisující síť \((x_n,\ n=0,1,\dots)\), y0 je počáteční podmínka \((y_{0r},\ r=1,\dots,R)\) a fcn(y,x) je funkce vyčíslující vektor pravých stran pro x = \(x\) a y = \((y_r,\ r=1,\dots,R)\). Cestou argumentu fcn lze dodat i funkci vracející Jacobiovu matici. Nastavení lsode (např. volba numerické metody) se děje voláním funkce lsode_options.

Příklad: 2D šikmý vrh \((x(t),z(t))\) bodu o hmotnosti \(M\) v gravitačním poli \(F_g=(0,-Mg)\) při odporové síle \(F_o=(-Kvv_x,-Kvv_z)\). Podle rovnic z ukázky v úvodu kapitoly volíme vektor neznámých v pořadí y = \((x,z,v_x,v_z)\).

Soubor oderhs.m s funkcí vyčíslující vektor pravých stran:

function result = oderhs(y,t);                   % funkce pro pravé strany (right-hand sides)
global M g K;                                    % data tekoucí mimo rozhraní
v=sqrt(y(3)^2+y(4)^2);                           % velikost rychlosti
result=[y(3); y(4); -K/M*v*y(3); -g-K/M*v*y(4)]; % (vx, vz, -K/M.v.vx, -K/M.v.vz)
end

Příkazy pro výpočet a vykreslení balistické křivky:

global M g K;
M=1; g=9.81; K=0.5; tvec=0:.1:2.5; y0=[0 0 7 7]; % parametry, sit, pocatecni podminka
y=lsode(@oderhs,y0,tvec);                        % Livermore Solver for ODEs
plot(y(:,1),y(:,2))

Totožnou úlohu řešíme i na stránkách cvičení (metoda, balistické křivky), kde je zpracován obdobný program v Pascalu a jsou vržena dvě reálná tělesa.