01NUM:Kapitola3
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 parabolického typu} Buďte $\Omega\subset\mathbbm R^n$ omezená oblast, jejíž hranicí je...)
[ 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 parabolického typu} Buďte $\Omega\subset\mathbbm R^n$ omezená oblast, jejíž hranicí je nadplocha $\Gamma$ po částech třídy $\mathcal C^1$, $L$ je eliptický parciální diferenciální operátor na $\Omega$, $T>0$. Zabývejme se smíšenou úlohou pro lineární parabolickou parciální diferenciální rovnici \[ \frac{\partial y}{\partial t}-Ly=f\quad\text{na }(0,\,T)\times\Omega \] společně s okrajovou podmínkou \[ \alpha(x)\frac{\partial y}{\partial\vec n}+\beta(x)y=\gamma(x)\quad\text{na }(0,\,T)\times\Gamma \] a počáteční podmínkou \[ y(0,\,x)=y_0(x)\quad\text{na }\Omega. \] \subsection{Metoda sítí} Omezíme se na následující jednorozměrný případ: \begin{subequations} \label{jednoduchaparabolickauloha} \begin{equation} \label{parabolickarce} \frac{\partial y}{\partial t}-D\frac{\partial^2y}{\partial x^2}=f(t,\,x)\quad\text{na }(0,\,T)\times(a,\,b), \quad D > 0 \end{equation} \begin{equation} \label{parabolickeop} y(t,\,a)=\gamma_1,\quad y(t,\,b)=\gamma_2\quad\text{na }(0,\,T), \end{equation} \begin{equation} y(0,\,x)=y_0(x)\quad\text{na }(a,\,b). \end{equation} \end{subequations} Buď $N_T\in\mathbbm N$ a označme $\tau=\frac T{N_T}$ časový krok. Na interval $(a,\,b)$ položíme prostorovou síť $\overline\omega_h$. Pro $k\in\{0,\,\hdots,\,N_T\}$, $j\in\{0,\,\hdots,\,m\}$ klademe \[ y^k_j:=y(k\tau,\,jh). \] Při řešení rovnic uvedeného typu vycházíme z formální analogie s obyčejnou diferenciální rovnicí: Místo toho, abychom na řešení pohlíželi jako na funkci $(t,\,x)\mapsto y(t,\,x)$, považujeme je za zobrazení $t\mapsto y(t,\,\cdot\,)$, které každému času $t\in(0,\,T)$ přiřazuje funkci prostorové proměnné. Síťovou funkci, která aproximuje řešení v~čase $t=k\tau$, resp. $t=(k+1)\tau$, budeme značit $u$, resp. $\widehat u$. \subsubsection{Explicitní schéma} Rovnici (\ref{parabolickarce}) nahradíme diferenční rovnicí \[ \frac{\widehat u-u}\tau-Du_{\bar xx}=f\quad\text{na }\omega_h. \] Odtud \[ \widehat u=u+\tau Du_{\bar xx}+\tau f, \] resp. po uzlech \[ u^{k+1}_j=u^k_j+\frac\tau{h^2}D\left(u^k_{j+1}-2u^k_j+u^k_{j-1}\right)+\tau f^k_j. \] Maticově můžeme získanou úlohu zapsat ve tvaru \begin{equation} \label{explicitnischema1} \widehat u=\A u+\tau f, \end{equation} kde \[ \A= \begin{pmatrix} 1-D\frac{2\tau}{h^2} & D\frac\tau{h^2} & 0 & 0 & \hdots\\ D\frac\tau{h^2} & 1-D\frac{2\tau}{h^2} & D\frac\tau{h^2} & 0 & \hdots\\ & \ddots & \ddots & \ddots \end{pmatrix}. \] Opakovanou aplikací (\ref{explicitnischema1}) dostaneme \[ u^k=\A u^{k-1}+\tau f^{k-1}=\A^2u^{k-2}+\tau\A f^{k-2}+\tau f^{k-1}=\hdots=\A^ku^0+\hdots. \] Protože si nepřejeme, aby malá změna $u^0$ mohla způsobit velkou změnu $u^k$, chceme, aby $\sigma(\A)\subset(-1,\,1)$. Matice $\A$ je 3-diagonální a má vlastní čísla \[ \lambda_i=1-4\frac\tau{h^2}D\sin^2\frac{i\pi}{2m},\quad i\in\widehat{m-1}. \] \begin{remark} Odpovídající vlastní vektory jsou \[ \left\{\sin\frac{i\pi j}m\right\}_{j=1}^{m-1},\quad i\in\widehat{m-1}. \] \end{remark} Požadujeme tedy, aby platilo \[ -1<1-4\frac\tau{h^2}D\sin^2\frac{i\pi}{2m}<1,\quad i\in\widehat{m-1}, \] tj. \[ \frac\tau{h^2}D\sin^2\frac{i\pi}{2m}<\frac12,\quad i\in\widehat{m-1}. \] Toho dosáhneme volbou \[ D\frac\tau{h^2}<\frac12. \] Explicitní schéma je proto {\em podmíněně stabilní}. \subsubsection{Implicitní schéma} Rovnici (\ref{parabolickarce}) nahradíme diferenční rovnicí \[ \frac{\widehat u-u}\tau-D\widehat u_{\bar xx}=\widehat f\quad\text{na }\omega_h. \] Odtud \[ \widehat u-\tau D\widehat u_{\bar xx}=u+\tau\widehat f, \] resp. po uzlech \[ u^{k+1}_j-D\frac\tau{h^2}\left(u^{k+1}_{j+1}-2u^{k+1}_j+u^{k+1}_{j-1} \right)=u^k_j+\tau f^{k+1}_j. \] Maticově můžeme získanou úlohu zapsat ve tvaru \begin{equation} \label{implicitnischema1} \A \widehat u=u+\tau\widehat f, \end{equation} kde \[ \A= \begin{pmatrix} 1+D\frac{2\tau}{h^2} & -D\frac\tau{h^2} & 0 & 0 & \hdots\\ -D\frac\tau{h^2} & 1+D\frac{2\tau}{h^2} & -D\frac\tau{h^2} & 0 & \hdots\\ & \ddots & \ddots & \ddots \end{pmatrix}. \] Opakovanou aplikací (\ref{implicitnischema1}) dostaneme \[ u^k=\A^{-1}(u^{k-1}+\tau f^k)=(\A^{-1})^2u^{k-2}+\tau\A^{-1}f^k+\tau(\A^{-1})^2f^{k-1}=\hdots=(\A^{-1} )^ku^0+\hdots. \] \begin{tvrz} Implicitní schéma je {\em nepodmíněně stabilní}. \begin{proof} Později dokážeme, že $\sigma(\A^{-1})\subset(-1,\,1)$. \end{proof} \end{tvrz} \subsubsection{Chyba aproximace diferenciálního operátoru $\frac\partial{\partial t}-L$ pro explicitní a implicitní schéma} Diferenciální operátor $y\mapsto\frac{\partial y}{\partial t}$ aproximujeme dopřednou diferencí $u_t$ v explicitním schématu a zpětnou diferencí $u_{\bar t}$ v implicitním schématu. Chyba je v obou případech řádu $O(\tau)$. Diferenciální operátor $y\mapsto-y''$ aproximujeme v obou schématech diferencí $u_{\bar xx}$. Chyba je řádu $O(h^2)$. Celková chyba aproximace diferenciálního operátoru $\frac\partial{\partial t}-L$ je tudíž v obou případech řádu $O(\tau+h^2)$. \subsubsection{Crankovo-Nicolsonovo schéma} Rovnici (\ref{parabolickarce}) nahradíme diferenční rovnicí \begin{equation} \label{ceesschema1} \frac{\widehat u-u}\tau-\frac D2\left(\widehat u_{\bar xx}+u_{\bar xx}\right)=\frac12\left(f+\widehat f\right)\quad\text{na }\omega_h. \end{equation} Odtud \[ \widehat u-\tau\frac D2\widehat u_{\bar xx}=u+\frac{\tau D}2u_{\bar xx}+\frac\tau2\left(\widehat f+f\right). \] \begin{remark} Schéma nahrazuje rovnici (\ref{parabolickarce}) v čase $\left(k+\frac12\right)\tau$; \begin{itemize} \item výraz $\frac1\tau\left(u^{k+1}-u^k\right)$ nahrazuje $\frac{\partial y}{\partial t}\left((k+\frac12)\tau,\,\cdot\,\right)$ s přesností řádu $O(\tau^2)$ (jde o centrální diferenci), \item výraz $\frac12\left(u^{k+1}_{\bar xx}+u^k_{\bar xx}\right)$ nahrazuje $\frac{\partial^2y}{\partial x^2}\left((k+\frac12)\tau,\,\cdot\,\right)$ s přesností řádu $O(\tau^2+h^2)$, \item výraz $\frac12\left(f^{k+1}+f^k\right)$ nahrazuje $f\left((k+\frac12)\tau,\,\cdot\,\right)$ s přesností řádu $O(\tau^2)$. \end{itemize} \end{remark} \begin{tvrz} Crankovo-Nicolsonovo schéma je nepodmíněně stabilní. \begin{proof} Položme \[\widehat u_{\overline t}=\frac{\widehat u-u}\tau,\quad\phi=\frac12\left(\widehat f+f\right).\] Pak můžeme (\ref{ceesschema1}) přepsat jako \[ \widehat u_{\overline t}-\phi=\frac D2\left(\widehat u_{\bar xx}+u_{\bar xx}\right). \] Tuto rovnici skalárně vynásobme výrazem $2\tau\widehat u_{\overline t}$. Využijeme-li Greenovu formuli, dostaneme díky okrajovým podmínkám (\ref{parabolickeop}) \[ \begin{split} 2\tau\nor{\widehat u_{\overline t}}_h^2-2\tau\left(\widehat u_{\overline t},\,\phi\right)_h&=\tau D\left(\widehat u_{\overline t},\,\widehat u_{\bar xx}+u_{\bar xx}\right)_h=D\left(\widehat u-u,\,(\widehat u+u)_{\bar xx}\right)_h=\\ =&-D\left((\widehat u-u)_{\bar x},\,(\widehat u+u)_{\bar x}\right]=-D\rnor{\widehat u_{\bar x}}^2+D\rnor{u_{\bar x}}^2. \end{split} \] Pomocí Schwarzovy a Youngovy nerovnosti dostaneme \[ 2\tau\nor{\widehat u_{\overline t}}_h^2+D\rnor{\widehat u_{\bar x}}^2-D\rnor{u_{\bar x}}^2=2\tau\left(\widehat u_{\overline t},\,\phi\right)_h\le2\tau\nor{\widehat u_{\overline t}}_h\nor\phi_h\le\tau\nor{\widehat u_{\overline t}}_h^2+\tau\nor\phi_h^2 \] neboli \[ \tau\nor{\widehat u_{\overline t}}_h^2+D\rnor{\widehat u_{\bar x}}^2-D\rnor{u_{\bar x}}^2\le\tau\nor\phi_h^2. \] Tím spíš ovšem platí \[ D\rnor{\widehat u_{\bar x}}^2-D\rnor{u_{\bar x}}^2\le\tau\nor\phi_h^2. \] Odtud \[ \begin{split} D\rnor{u^{k+1}_{\bar x}}^2&\le D\rnor{u^k_{\bar x}}^2+\tau\nor{\phi^k}_h^2\le D\rnor{u^{k-1}_{\bar x}}^2+\tau\rnor{\phi^{k-1}}_h^2+\tau\nor{\phi^k}_h^2\le\hdots\le\\ &\le D\rnor{u^0_{\bar x}}^2+\sum_{j=0}^k\tau\nor{\phi^j}_h^2. \end{split} \] Původní úlohu (\ref{jednoduchaparabolickauloha}) nahrazujeme diferenční úlohou \begin{subequations} \label{nahradajpu} \begin{align} \frac{u^{k+1}-u^k}\tau&=\frac12D\left(u^{k+1}_{\bar xx}+u^k_{\bar xx}\right)+\frac12\left(f^{k+1}-f^k\right),&k&=1,\,\hdots,\,N_T-1,\\ u^k_0&=\gamma_1,\quad u^k_m=\gamma_2,&k&=1,\,\hdots,\,N_T,\\ u^0_j&=y_{0j},&j&=0,\,\hdots,\,m. \end{align} \end{subequations} Odečtěme spojitou úlohu pro hladinu $k+\frac12$ a tuto diskrétní úlohu. Chyba aproximace diferenciálního operátoru $\frac\partial{\partial t}-D\frac{\partial^2}{\partial x^2}$ je \[ \begin{split} \Psi_{h,\tau}&=\mathcal P_{h,\tau} \left(\frac{\partial y}{\partial t}-D\frac{\partial^2 y}{\partial x^2}\right)- \frac{\left(\mathcal P_{h,\tau}y\right)^{k+1}-\left(\mathcal P_{h,\tau}y\right)^k}\tau+\\ &\quad+\frac12 D\left((\mathcal P_{h,\tau}y)_{\bar xx}^{k+1}+ (\mathcal P_{h,\tau}y)_{\bar xx}^{k} \right)=O(h^2+\tau^2), \end{split} \] kde $(\mathcal P_{h,\tau}y)_j^k=y(k\tau,\,a+jh)$, $j=0,\dots,m$, $k=0,\dots,N_T$. Položme $z=u-\mathcal P_{h,\tau}y$. Potom $z$ splňuje soustavu rovnic \[ \begin{split} \frac{\widehat z-z}{\tau}&=\frac12D(\widehat z_{\bar xx}+z_{\bar xx})+\Psi_{h,\tau},\\ z_0^k=z_m^k&=0,\quad k=1,\,\dots,\,N_T,\\ z_j^0&=0,\quad j=0,\,\dots,\,m. \end{split} \] Odtud a z předchozího vyplývá, že platí \[ D\rnor{z^{k+1}_{\bar x}}^2\le D\rnor{z^0_{\bar x}}^2+\sum_{j=0}^k\tau\nor{\Psi^j}_h^2=\sum_{j=0}^k\tau\nor{\Psi^j}_h^2. \] Ze Sobolevových nerovností pak plyne \[ \nor{z^{k+1}}_{h,0}\le\frac{\sqrt{b-a}}2\rnor{z^{k+1}_{\bar x}}\le\frac{\sqrt{b-a}}{2D}\sqrt{\sum_{j=0}^k\tau\nor{\Psi^j}_h^2}=O(\tau^2+h^2) , \] %\renewcommand{\qed}{} tj. stabilita a konvergence. \end{proof} \end{tvrz} \subsection{Metoda přímek} V metodě sítí jsme vycházeli z toho, že řešení úlohy (\ref{jednoduchaparabolickauloha}) je zobrazení, které každému časovému okamžiku přiřadí funkci prostorové proměnné. Nabízí se též opačný pohled na věc, tj. sledovat časový vývoj ve zvoleném bodě $x\in\langle a,\,b\rangle$ jako funkci času. Postupujeme tak, že na interval $\langle a,\,b\rangle$ položíme síť a úlohu (\ref{jednoduchaparabolickauloha}) převedeme na počáteční úlohu pro soustavu obyčejných diferenciálních rovnic \begin{align*} \frac{\text du_j}{\text dt}&=D\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}+f_j,&j\in\widehat{m-1},\\ u_0(t)&=\gamma_1,\quad u_m(t)=\gamma_2,&t\in(0,\,T),\\ u_j(0)&=y_{0j},&j\in\widehat{m-1}, \end{align*} stručně zapsáno \begin{align*} \frac{\text d\tucne u}{\text dt}&=D\tucne u_{\bar xx}+\tucne f,\\ \tucne u(0)&=\tucne y_0. \end{align*} Tuto soustavu řešíme např. Rungovou-Kuttovou metodou.