Zdrojový kód
%\wikiskriptum{01NUM}
\section{Numerické řešení okrajových úloh pro PDE eliptického typu}
Buď $\Omega\subset\mathbbm R^n$ omezená oblast, jejíž hranicí je nadplocha
$\Gamma$ po částech třídy $\mathcal C^1$. Zabývejme se lineární parciální
diferenciální rovnicí 2. řádu
\begin{equation}
-\diverge(\matice A(x)\nabla y)+q(x)y=f(x)\quad\text{v }\Omega
\end{equation}
společně s okrajovou podmínkou
\begin{equation}
\alpha(x)\frac{\partial y}{\partial\vec n}+\beta(x)y=\gamma(x)\quad\text{na }\Gamma.
\end{equation}
Symbol $\frac{\partial y}{\partial\vec n}$ značí derivaci ve směru tzv. {\em
konormály} $\vec n=\matice A^{\text{\sf T}}\vec\nu$, kde $\matice A^{\text{\sf
T}}$ značí transponovanou matici $\matice A$ a $\vec\nu$ je směrový vektor
vnější normály k nadploše $\Gamma$.
Předpokládáme, že platí $q\ge0$ a že existuje $p_0>0$ tak, že
\[
\left(\matice A\vec\xi,\,\vec\xi\right)\ge p_0\nor{\vec\xi}^2\quad\text{na}\Omega,\quad\forall\vec\xi\in\mathbbm R^n.
\]
\begin{remark}
Jsou-li $\alpha,\,\beta\ge0$ a navíc $\alpha+\beta>0$ a jsou-li funkce
$\A,\,q,\,f,\,\alpha,\,\beta,\,\gamma$ dost hladké, je úloha jednoznačně
řešitelná.
\end{remark}
\subsection{Metoda sítí}
Omezíme se na obdélníkovou oblast v $\mathbbm R^2$, tj. bez újmy na obecnosti
oblast tvaru $\Omega=(0,\,L_1)\times(0,\,L_2)$. Buďte $m_1,\,m_2\in\mathbbm N$ a
položme
\[
h_1=\frac{L_1}{m_1},\quad h_2=\frac{L_2}{m_2}.
\]
Na $\overline\Omega$ položíme síť
\[
\overline\omega_h=\left\{[ih_j,\,jh_2]\,|\,i=0,\,\hdots,\,m_1;\,j=0,\,\hdots,\,m_2\right\}.
\]
Množina vnitřních bodů sítě je
\[
\omega_h=\left\{[ih_j,\,jh_2]\,|\,i=1,\,\hdots,\,m_1-1;\,j=1,\,\hdots,\,m_2-1\right\}
\]
a množina hraničních uzlů
\[
\gamma_h=\overline\omega_h-\omega_h.
\]
Zabývejme se nyní úlohou
\begin{subequations}
\label{eliptickaulohapdr}
\begin{align}
-\frac\partial{\partial x_1}\left(p(x)\frac{\partial y}{\partial x_1}\right)
-\frac\partial{\partial x_2}\left(p(x)\frac{\partial y}{\partial x_2}\right)
+q(x)y&=f(x)\quad\text{na }\Omega,\\
y&=\gamma_0\quad\text{na }\Gamma,
\end{align}
\end{subequations}
kde $\gamma_0\in\mathbbm R$.
\begin{remark}
Hodnoty funkcí v síťových uzlech značíme $y_{ij}=y(ih_1,\,jh_2)$,
$i\in\widehat{m_1}_0$, $j\in\widehat{m_2}_0$.
Parciální derivace funkcí aproximujeme pomocí diferencí
\[
\left(\frac{\partial y}{\partial x_1}\right)_{ij}\approx y_{\bar x_1}=\frac{y_{ij}-y_{i-1j}}{h_1},\quad
\left(\frac{\partial y}{\partial x_2}\right)_{ij}\approx y_{\bar x_2}=\frac{y_{ij}-y_{ij-1}}{h_2}.
\]
Přesnost je v obou případech řádu $O(h_1+h_2)$.
\end{remark}
Úlohu (\ref{eliptickaulohapdr}) nahradíme diferenčním schématem
\begin{subequations}
\label{eliptickadiferencniuloha}
\begin{align}
\label{eliptickediferencischema}
-(pu_{\bar x_1})_{x_1}-(pu_{\bar x_2})_{x_2}+qu&=f\quad\text{na }\omega_h,\\
u&=\gamma_0\quad\text{na }\gamma_h.
\end{align}
\end{subequations}
\begin{remark}[5bodové schéma]
V celé této poznámce bude index $i$, resp. $j$ probíhat množinu
$\widehat{m_1-1}$, resp. $\widehat{m_2-1}$. Podle předchozí poznámky můžeme
rozepsat rovnici (\ref{eliptickediferencischema}) jako
\begin{multline*}
-\frac1{h_1}\left(p_{i+1j}\frac{u_{i+1j}-u_{ij}}{h_1}-p_{ij}\frac{u_{ij}-u_{i-1j}}{h_1}\right)-\\
-\frac1{h_2}\left(p_{ij+1}\frac{u_{ij+1}-u_{ij}}{h_2}-p_{ij}\frac{u_{ij}-u_{ij-1}}{h_2}\right)
+q_{ij}u_{ij}=f_{ij}
\end{multline*}
neboli
\[
A_{ij}u_{i-1j}+B_{ij}u_{ij-1}+C_{ij}u_{i+1j}+D_{ij}u_{ij+1}+E_{ij}u_{ij}=F_{ij},\]
kde
\[
\begin{split}
A_{ij}&=-\frac{p_{ij}}{h_1^2},\quad B_{ij}=-\frac{p_{ij}}{h_2^2},\quad C_{ij}=-\frac{p_{i+1j}}{h_1^2},\quad D_{ij}=-\frac{p_{ij+1}}{h_2^2},\\
E_{ij}&=\frac{p_{i+1j}}{h_1^2}+\frac{p_{ij}}{h_1^2}+\frac{p_{ij}}{h_2^2}+\frac{p_{ij+1}}{h_2^2}+q_{ij},\quad F_{ij}=f_{ij}.
\end{split}
\]
Tuto soustavu lze řešit např. zobecněnou faktorizací.
\end{remark}
\subsection{Konvergence, odhad chyby}
Restrikci funkce $y:\overline{\Omega}\rightarrow\mathbbm R$ na síť $\omega_h$
budeme opět značit $\mathcal P_h y$. Dále zavádíme chybu aproximace
diferenciálního operátoru $L$ diferenčním operátorem $L_h$ jako síťovou funkci
$\Psi_h:\overline\omega_h\rightarrow\mathbbm R$, $\Psi_h=\mathcal
P_h(Ly)-L_h(\mathcal P_h y)$.
\begin{define}
Buďte $u,\,v:\overline\omega_h\rightarrow\mathbbm R$. Potom definujeme skalární
součiny
\[
\begin{split}
(u,v)_h&=\sum_{i=1}^{m_1-1}\sum_{j=1}^{m_2-1}h_1h_2u_{ij}v_{ij},\\
(u,v\rfloor&=\sum_{i=1}^{m_1}\sum_{j=1}^{m_2-1}h_1h_2u_{ij}v_{ij},\\
(u,v\rceil&=\sum_{i=1}^{m_1-1}\sum_{j=1}^{m_2}h_1h_2u_{ij}v_{ij}.
\end{split}
\]
\end{define}
\begin{define}
Buď $u:\overline\omega_h\rightarrow\mathbbm R$. Potom definujeme normu
\[\nor{u}_h=\sqrt{(u,\,u)_h}.\]
\end{define}
\begin{define}
Buďte $U,\,V:\overline\omega_h\rightarrow\mathbbm R^2$, $U=[U^1,\,U^2]$,
$V=[V^1,\,V^2]$. Potom klademe
\[(U,V]=(U^1,V^1\rfloor+(U^2,V^2\rceil.\]
\end{define}
\begin{define}
Buď $U:\overline\omega_h\rightarrow\mathbbm R^2$, $U=[U^1,\,U^2]$. Potom
definujeme normu
\[\rnor{U}=\sqrt{(U,\,U]}.\]
\end{define}
\begin{define}
Buď $u:\overline\omega_h\rightarrow\mathbbm R$. Potom definujeme síťový a zpětný
gradient
\[
\nabla_hu=[u_{x_1},\,u_{x_2}],\quad\onabla_hu=[u_{\bar x_1},\,u_{\bar x_2}],
\]
a síťový laplacián
\[
\Delta_hu=u_{\bar x_1x_1}+u_{\bar x_2x_2}.
\]
\end{define}
\begin{define}
Buď $U:\overline\omega_h\rightarrow\mathbbm R^2$, $U=[U^1,\,U^2]$. Potom
definujeme síťovou divergenci
\[\diverge_h U=(U^1)_{x_1}+(U^2)_{x_2}.\]
\end{define}
\begin{tvrz}[Greenova formule]
Buďte $u,\,v:\overline\omega_h\rightarrow\mathbbm R$ a nechť $u=v=0$ na
$\gamma_h$. Pak platí
\[\left(\diverge_h(p\onabla_hu),\,v\right)_h=-\left(p\onabla_hu,\,
\onabla_hv\right].\]
\begin{proof}
K dispozici máme jednorozměrnou Greenovu formuli: Jsou-li $u,\,v$ funkce
definované na jednorozměrné síti a splňující $u_0=u_m=v_0=v_m=0$, pak
\[(v,(pu_{\bar x})_x)_h=-(pu_{\bar x},v_{\bar x}].\]
Tento vztah můžeme použít pro jednorozměrná zúžení síťových funkcí $u$, $v$ ve
směru $x_1$. Dostaneme
\[\left(v_{\cdot j},{(pu_{\bar x_1})_{x_1}}_{\cdot j}\right)_h=-((pu_{\bar x_1})_{\cdot j},\,{v_{\bar x_1}}_{\cdot j}],\quad j\in\widehat{m_2-1},\]
neboli
\[\sum_{i=1}^{m_1-1}h_1v_{ij}{(pu_{\bar x_1})_{x_1}}_{ij}=-\sum_{i=1}^{m_1}h_1p_{ij}{u_{\bar x_1}}_{ij}{v_{\bar x_1}}_{ij},\quad j\in\widehat{m_2-1}.\]
Vynásobením $h_2$ a vysčítáním přes $j$ dostáváme
\[(v,(pu_{\bar x_1})_{x_1})_h=-(v_{\bar x_1},\,pu_{\bar x_1}\rfloor.\]
Obdobně bychom odvodili vztah
\[(v,(pu_{\bar x_2})_{x_2})_h=-(v_{\bar x_2},\,pu_{\bar x_2}\rceil.\]
Sečtením posledních dvou rovností pak dostaneme
\[
\begin{split}
\left(\diverge_h(p\onabla_hu),\,v\right)_h&=\left((pu_{\bar x_1})_{x_1}+(pu_{\bar x_2})_{x_2},\,v\right)_h=\\
&=-(v_{\bar x_1},\,pu_{\bar x_1}\rfloor-(v_{\bar x_2},\,pu_{\bar x_2}\rceil=-\left(\onabla_h v,\,p\onabla_h u\right].\qed
\end{split}
\]
\renewcommand{\qed}{}
\end{proof}
\end{tvrz}
\begin{lemma}[Sobolev]
\label{sobolev2_2d}
Nechť $u:\overline\omega_h\rightarrow\mathbbm R$, $u=0$ na $\gamma_h$. Pak
\[
\nor{u}_h^2\le\frac18\max\{L_1^2,\,L_2^2\}\rnor{\onabla_hu}^2.
\]
\end{lemma}
\begin{proof}
Podle lemmatu \ref{sobolev2}, aplikovaného na jednorozměrná zúžení funkce $u$
ve směru osy $x_1$, resp. $x_2$ platí
\[\nor{u_{\cdot j}}_h\le\frac{L_1}{2}\rnor{(u_{\cdot j})_{\bar x_1}},\quad\forall j\in\widehat{m_2-1},\]
resp.
\[\nor{u_{i\cdot}}_h\le\frac{L_2}{2}\rnor{(v_{i\cdot})_{\bar x_2}},\quad\forall i\in\widehat{m_1-1}.\]
Po umocnění na druhou můžeme tyto odhady přepsat jako
\begin{gather}
\label{eq:sob2d_1}
\sum_{i=1}^{m_1-1}h_1\abs{u_{ij}}^2\le \frac{L_1^2}4\sum_{i=1}^{m_1}h_1\abs{{u_{\bar x_1}}_{ij}}^2, \quad\left|\cdot h_2,\sum_{j=1}^{m_2-1}\right.\\
\label{eq:sob2d_2}
\sum_{j=1}^{m_2-1}h_1\abs{u_{ij}}^2\le \frac{L_2^2}4\sum_{i=1}^{m_2}h_2\abs{{u_{\bar x_2}}_{ij}}^2. \quad\left|\cdot h_1,\sum_{i=1}^{m_1-1}\right.
\end{gather}
Sečtením \eqref{eq:sob2d_1} a \eqref{eq:sob2d_2} dostaneme
\[
\begin{split}
\nor{u}_h^2&=\sum_{i=1}^{m_1-1}\sum_{j=1}^{m_2-1}h_1h_2\abs{u_{ij}}^2\le\\
&\le\frac{L_1^2}8\sum_{i=1}^{m_1}\sum_{j=1}^{m_2-1}h_1h_2 \abs{{u_{\bar x_1}}_{ij}}^2+ \frac{L_2^2}8\sum_{i=1}^{m_1-1}\sum_{j=1}^{m_2}h_1h_2 \abs{{u_{\bar x_2}}_{ij}}^2\le\\
&\le\frac18\max\left\{L_1^2,\,L_2^2\right\}\left(\sum_{i=1}^{m_1}\sum_{j=1}^{m_2-1}h_1h_2\abs{{u_{\bar x_1}}_{ij}}^2+\sum_{i=1}^{m_1-1}\sum_{j=1}^{m_2}h_1h_2\abs{{u_{\bar x_2}}_{ij}}^2\right)=\\
&=\frac18\max\left\{L_1^2,\,L_2^2\right\}\rnor{\onabla_h u}^2.\qed
\end{split}
\]
\renewcommand{\qed}{}
\end{proof}
\subsubsection{Metoda energetických nerovností}
Úlohu (\ref{eliptickaulohapdr}) zúžíme na síť a odečteme ji od
(\ref{eliptickadiferencniuloha}). Položíme-li
\[
L:y\mapsto-\frac\partial{\partial x_1}\left(p(x)\frac{\partial y}{\partial x_1}\right)-\frac\partial{\partial x_2}\left(p(x)\frac{\partial y}{\partial x_2}\right)+q(x)y,\]
$L_h:u\mapsto-(pu_{\bar x_1})_{x_1}-(pu_{\bar x_2})_{x_2}+qu$, dostaneme
\begin{subequations}
\label{energetickyprepiseliptickeulohymain}
\begin{align}
\label{energetickyprepiseliptickeulohyrce}
L_hu-\mathcal P_h(Ly)&=0\quad\text{na }\omega_h,\\
u-\mathcal P_hy&=0\quad\text{na }\gamma_h.
\end{align}
\end{subequations}
Chyba aproximace $\Psi_h=\mathcal P_h(Ly)-L_h(\mathcal P_hy)$ je řádu
$O(h_1+h_2)$ a s její pomocí můžeme (\ref{energetickyprepiseliptickeulohyrce})
přepsat jako
\[
L_hu-L_h(\mathcal P_hy)=\Psi_h.
\]
Položíme-li $z=u-\mathcal P_hy$, můžeme tudíž
(\ref{energetickyprepiseliptickeulohymain}) psát ve tvaru
\begin{align}
L_hz&=\Psi_h\quad\text{na }\omega_h,\\
z&=0\quad\text{na }\gamma_h.
\end{align}
První z obou vztahů skalárně vynásobíme $z$. S využitím Greenovy formule pak
dostaneme
\[
\begin{split}
(\Psi_h,\,z)_h&=(L_h z,\,z)_h=\left(-(pz_{\bar x_1})_{x_1}-(pz_{\bar x_2})_{x_2}+(qz,\,z\right)_h=\\
&=-\left(\diverge_h(p\onabla_hz),\,z\right)_h+(qz,\,z)_h=\left(p\onabla_hz,\,\onabla_hz\right]+(qz,\,z)_h=\\
&=(pz_{\bar x_1},\,z_{\bar x_1}\rfloor+(pz_{\bar x_2},\,z_{\bar x_2}\rceil+(qz,\,z)_h.
\end{split}
\]
Podle základních předpokladů je $q\ge 0$, a proto $(qz,\,z)_h\ge 0$. Dále je
$p\ge p_0>0$, takže
\[\left(p\onabla_hz,\,\onabla_hz\right]\ge p_0\rnor{\onabla_hz}^2.\]
Celkově jsme dokázali, že
\[(\Psi_h,\,z)_h\ge p_0\rnor{\onabla_hz}^2.\]
Podle lemmatu \ref{sobolev2_2d} a Schwarzovy nerovnosti je
\[
\nor{z}_h^2\le c\rnor{\onabla
z}^2\le\frac{c}{p_0}(\Psi_h,\,z)_h\le\frac{c}{p_0}\nor{\Psi_h}_h\nor{z}_h
\]
a odtud
\[\nor{z}_h\le\frac{c}{p_0}\nor{\Psi_h}_h.\]
Diferenční schéma (\ref{eliptickadiferencniuloha}) je tedy korektní (mj. také
stabilní) a $\nor{z}_h$ se chová stejně jako $\nor{\Psi_h}_h$, tj. jako
$O(h_1+h_2)$.
\begin{remark}
V případě, že $p\equiv\text{konst.}$, je dokonce $\Psi_h=O(h_1^2+h_2^2)$.
\end{remark}
