01NUM:Kapitola2
Z WikiSkripta FJFI ČVUT v Praze
Verze z 1. 8. 2010, 01:49, kterou vytvořil Admin (diskuse | příspěvky) (Založena nová stránka: %\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 nad...)
[ znovu generovat, | výstup z překladu ] | Kompletní WikiSkriptum včetně všech podkapitol. | |
PDF Této kapitoly | [ znovu generovat, | výstup z překladu ] | Přeložení pouze této kaptioly. |
ZIP | Kompletní zdrojový kód včetně obrázků. |
Součásti dokumentu 01NUM
součást | akce | popis | poslední editace | soubor | |||
---|---|---|---|---|---|---|---|
Hlavní dokument | editovat | Hlavní stránka dokumentu 01NUM | Johndavi | 3. 6. 2019 | 16:06 | ||
Řídící stránka | editovat | Definiční stránka dokumentu a vložených obrázků | Admin | 7. 9. 2015 | 14:49 | ||
Header | editovat | Hlavičkový soubor | Johndavi | 1. 6. 2019 | 18:30 | header.tex | |
Kapitola1 | editovat | Numerické řešení okrajových úloh pro ODE | Kunzmart | 2. 9. 2021 | 00:00 | kapitola1.tex | |
Kapitola2 | editovat | Numerické řešení okrajových úloh pro PDE eliptického typu | Kunzmart | 1. 9. 2021 | 23:28 | kapitola2.tex | |
Kapitola3 | editovat | Numerické řešení okrajových úloh pro PDE parabolického typu | Admin | 1. 8. 2010 | 01:49 | kapitola3.tex | |
Kapitola4 | editovat | Numerické řešení okrajových úloh pro PDE 1. řádu | Johndavi | 3. 6. 2019 | 19:59 | kapitola4.tex | |
Kapitola5 | editovat | Numerické řešení počátečních úloh pro ODE | Johndavi | 30. 4. 2019 | 00:47 | kapitola5.tex |
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}