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

(rozdíl) ← Starší verze | zobrazit aktuální verzi (rozdíl) | Novější verze → (rozdíl)
Přejít na: navigace, hledání
PDF [ 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.
ZIPKompletní zdrojový kód včetně obrázků.

Součásti dokumentu 01NUM

součástakcepopisposlední editacesoubor
Hlavní dokument editovatHlavní stránka dokumentu 01NUMJohndavi 3. 6. 201916:06
Řídící stránka editovatDefiniční stránka dokumentu a vložených obrázkůAdmin 7. 9. 201514:49
Header editovatHlavičkový souborJohndavi 1. 6. 201918:30 header.tex
Kapitola1 editovatNumerické řešení okrajových úloh pro ODEKunzmart 2. 9. 202100:00 kapitola1.tex
Kapitola2 editovatNumerické řešení okrajových úloh pro PDE eliptického typuKunzmart 1. 9. 202123:28 kapitola2.tex
Kapitola3 editovatNumerické řešení okrajových úloh pro PDE parabolického typuAdmin 1. 8. 201001:49 kapitola3.tex
Kapitola4 editovatNumerické řešení okrajových úloh pro PDE 1. řáduJohndavi 3. 6. 201919:59 kapitola4.tex
Kapitola5 editovatNumerické řešení počátečních úloh pro ODE Johndavi 30. 4. 201900: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.

Citováno z „https://wikiskripta.fjfi.cvut.cz/wiki/index.php?title=01NUM:Kapitola3&oldid=3127“