Programování pro fyziky 2023/24 – Numerické recepty č. 5

Při výkladu si pomáháme odkazy do knih Numerical Recipes in C++, Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR) a také Numerical Computing with MATLAB, které jsou na webu volně k dispozici. Výklad provázejí ukázky skriptů pro Python s balíčky NumPy a SciPy a také pro Matlab a GNU Octave. V těchto prostředích je k demonstraci numerických metod krátká cesta.

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 (Newtonův druhý pohybový zákon): 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 pro naši úlohu 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 sobě 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\).

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

# Python: přesnost a stabilita explicitní a implicitní Eulerovy metody
import numpy as np
import matplotlib.pyplot as plt
nmax=5; ye=np.empty(nmax+1); yi=np.empty(nmax+1)
c=1; x0=0; y0=1; ye[0]=y0; yi[0]=y0; xmax=10.4; h=xmax/nmax; x=np.linspace(x0,xmax,nmax+1)
for n in range(nmax): ye[n+1]=(1-c*h)*ye[n]; yi[n+1]=yi[n]/(1+c*h)
plt.plot(x,np.exp(-c*x),'g',x,ye,'r',x,yi,'b')
plt.show()
% Octave: přesnost a stabilita explicitní a implicitní Eulerovy metody
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 tečen), 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 je tedy explicitní aproximací implicitní Eulerovy metody.

Rungovy-Kuttovy metody

(NR kap 16.1) Obecný vzorec (zde: explicitních) Rungových-Kuttových 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) vděčí za oblibu svým zázračně jednoduchým koeficientům, z nichž několik je dokonce nulových:

\[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 (také na Wikipedii).

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})\).

Python: SciPy a subbalíček integrate. Níže používáme řešič odeint. 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 neboli balistická křivka \((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)\). Vektor neznámých volíme v pořadí y = \((x,z,v_x,v_z)\), v souladu s pořadím rovnic v ukázce z úvodu kapitoly, tedy

\[\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}\]
# Python: řešení soustavy obyčejných diferenciálních rovnic pro 2D šikmý vrh
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
M=1; g=9.81; K=0.5                                 # fyzikální parametry
tvec=np.linspace(0,2.5,26); y0=[0,0,7,7]           # síť, počáteční podmínka

def oderhs(y,t):                                   # funkce pro pravé strany (right-hand sides)
  (x,z,vx,vz)=y
  v=np.sqrt(vx**2+vz**2)                           # velikost rychlosti
  return [vx,vz,-K/M*v*vx,-g-K/M*v*vz]

y=sp.integrate.odeint(oderhs,y0,tvec)              # solver for ODEs
plt.plot(y[:,0],y[:,1])
plt.show()
% Octave: řešení soustavy obyčejných diferenciálních rovnic pro 2D šikmý vrh
0;

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

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))
pause