Numerická matematika 2 ⬩ Adamátorovy zápisky

Hlavní téma: numerické řešení okrajových a smíšených úloh pro diferenciální rovnice

Kdo odevzdá praktickou část do určitého termínu (7. 6. 2024), bude mít u zkoušky o jednu otázku méně (3 → 2)

Metoda střelby

Jako jednoduchou typovou úlohu vezmeme jednu rovnici druhého řádu s Dirichletovými okrajovými podmínkami:

y^{″} = f (x, y, y^{'}), x \in (a, b), y (a) = γ_{1}, y (b) = γ_{2}

Myšlenka je taková, že si zvolíme hodnotu derivace na jednom okraji a budeme úlohu řešit jako počáteční. Jakmile se dostaneme na konec neboli na druhý okraj, nejspíš zjistíme, že jsme se dostali jinam, než jsme chtěli. Tak zvolíme jinou počáteční derivaci a zkusíme znovu.

Formálněji: Budeme opakovaně řešit počáteční úlohu (pomocí Eulerovy nebo Runge-Kuttovy metody)

y^{″} = f (x, y, y^{'}), y (a) = γ_{1}, y^{'} (a) = α

Tím dostaneme nějaké řešení $y (x; α)$ a chceme najít $α$ takové, aby $y (b; α) = γ_{2}$ . To můžeme hledat například pomocí metody půlení intervalu nebo Newtonovy metody.

Půlení intervalu tedy bude fungovat tak, že najdeme ${\bar{α}}_{1}, {\bar{α}}_{2}$ tak, aby $y (b; {\bar{α}}_{1}) < γ_{2}$ a $y (b; {\bar{α}}_{2}) > γ_{2}$ . Nastavíme $α_{1} \leftarrow {\bar{α}}_{1}, α_{2} \leftarrow {\bar{α}}_{2}, α ≔ \frac{α_{1} + α_{2}}{2}$ . Je-li $y (b; α) < γ_{2}$ , do další iterace vezmeme $α_{1} \leftarrow α$ , jinak $α_{2} \leftarrow α$ . Končíme, jakmile vzdálenost bude menší než nějaká tolerance. Většinou se používá relativní tolerance $10^{- 6}$ .

Ale pomocí Newtonovy metody to často konverguje rychleji. Rovnici $y (b; α) = γ_{2}$ můžeme zapsat jako $F (α) = 0$ , kde $F (α) ≔ y (b; α) - γ_{2}$ . Potom stačí iterovat podle vzorečku $α_{k + 1} ≔ α_{k} - {(F^{'} (α_{k}))}^{- 1} F (α_{k})$ . Akorát si musíme napočítat derivaci:

F^{'} (α) = \partial_{α} (y (b; α) - γ_{2}) = \partial_{α} y (b; α)

Tuto derivaci určíme řešením úlohy ve variacích.

Nebo pokud se nám nechce hledat derivaci tímhle způsobem, můžeme ji zkusit aproximovat podle definice, ovšem nic nám nezaručuje, že to bude dobrá aproximace:

F^{'} (α_{k}) \approx \frac{F (α_{k}) - F (α_{k - 1})}{α_{k} - α_{k - 1}}

Newtonova metoda poskytuje posloupnost $(α_{k})$ konvergující k řešení, ale $| F (α_{k}) |$ nemusí konvergovat monotónně. Pokud nám to vadí, můžeme použít modifikaci, jejíž myšlenka je v tom, že se budeme pohybovat o menší kus, abychom nepřestřelili:

α_{k + 1} ≔ α_{k} - λ_{k} {(F^{'} (α_{k}))}^{- 1} F (α_{k}), λ_{k} \in (0, 1 ⟩

Samozřejmě je nutné správně zvolit $λ_{k}$ . Označme

Δ α_{k} ≔ - {(F^{'} (α_{k}))}^{- 1} F (α_{k}) (o kolik se posouváme při normální Newtonově metodě)

Φ (α) ≔ {‖ F (α) ‖}_{2}^{2} = \sum_{j = 1}^{n} F_{j} {(α)}^{2} (jakási chyba aproximace)

φ (λ) ≔ Φ (α_{k} + λ Δ α_{k}) (chyba, když se posuneme o $\lambda$-násobek normálního kroku)

Spočteme si

φ^{'} (λ) = 2 \sum_{j = 1}^{n} F_{j} (α_{k} + λ Δ α_{k}) \cdot \sum_{l = 1}^{n} \partial_{α_{l}} F_{j} (α_{k} + λ Δ α_{k}) \cdot {(Δ α_{k})}_{l} = 2 F {(α_{k} + λ Δ α_{k})}^{𝖳} \cdot (F^{'} (α_{k} + λ Δ α_{k}) Δ α_{k})

φ^{'} (0) = 2 F {(α_{k})}^{𝖳} \cdot (F^{'} (α_{k}) Δ α_{k}) = - 2 F {(α_{k})}^{𝖳} \cdot (F^{'} (α_{k}) {(F^{'} (α_{k}))}^{- 1} F (α_{k})) = - 2 Φ (α_{k}) < 0

Jelikož $φ \in 𝒞^{1}$ , jistě najdeme $λ_{0} \in ℝ^{+}$ takové, že $φ^{'} (λ) < 0$ pro $λ \in ⟨ 0, λ_{0})$ .

Algoritmus bude vypadat tak, že nejprve zvolíme $λ_{k} = 1$ . Pokud $Φ (α_{k} + λ_{k} Δ α_{k}) < Φ (α_{k})$ , spočteme si $α_{k + 1} ≔ α_{k} + λ_{k} Δ α_{k}$ a budeme normálně pokračovat. Pokud ne, pokusíme se najít vhodné $λ$ – například budeme opakovaně půlit, až bude nerovnost platit (což z teorie plyne, že nastane). Každopádné skončíme, až bude $Φ (α_{k}) < ε^{2}$ , kde $ε$ je zadaná tolerance.

Metoda střelby pro okrajovou podmínku nelineárního typu

Typová úloha:

y^{'} = f (x, y), r (y (a), y (b)) = 0, f : ℝ^{1 + n} \to ℝ^{n}, r : ℝ^{n + n} \to ℝ^{n}, n \geq 2

Parametry funkce $r$ budeme značit $ξ, η$ .

y_{1^{'}} = y^{2}

Budeme řešit podobnou počáteční úlohu, ale s tím, že $α$ je tentokrát vektor:

y^{'} = f (x, y), y (a) = α

Chceme najít správné $α$ pomocí okrajové podmínky. Budeme tedy řešit rovnici

r (α, y (b; α)) = 0

Tu budeme řešit pomocí Newtonovy metody, kde vezmeme $F (α) ≔ r (α, y (b; α))$ .

F^{'} {(α_{k})}_{i, j} = \partial_{α_{j}} F_{i} (α_{k}) = \partial_{ξ_{j}} r_{i} (α_{k}, y (b; α_{k})) + \sum_{l = 1}^{n} \partial_{η_{l}} r_{i} (α_{k}, y (b; α_{k})) \cdot \partial_{α_{j}} y_{l} (b; α_{k})

Ke spočtení tohoto potřebujeme řešit soustavu ve variacích:

\partial_{x} \partial_{α_{j}} y_{i} (x; α) = \sum_{l = 1}^{n} \partial_{y_{l}} f_{i} (x, y (x; α)) \cdot \partial_{α_{j}} y_{l} (x; α), \partial_{α_{j}} y_{i} (a; α) = δ_{i, j}

Metoda střelby pro soustavu lineárních rovnic

y^{'} = 𝐀 (x) y + b (x), 𝐔 y (a) + 𝐕 y (b) = c

𝐀 : ⟨ a, b ⟩ \to ℝ^{n \times n}, 𝒞^{0}

b : ⟨ a, b ⟩ \to ℝ^{n}, 𝒞^{0}

(

b

jako

b

eneš není schopný značit různé věci různými písmeny)

c \in ℝ^{n}, 𝐔, 𝐕 \in ℝ^{n \times n}, h (𝐔) + h (𝐕) = n

Opakovaně řešíme počáteční úlohu

y^{'} = 𝐀 (x) y + b (x), y (a) = α

Využijeme toho, že kdyby bylo jen $y^{'} = 𝐀 (x) y$ , dokážeme soustavu řešit analyticky pomocí matice fundamentálního systému:

Φ (x) ≔ [\begin{matrix} v_{1}^{1} (x) & \dots & v_{n}^{1} (x) \\ ⋮ & ⋱ & ⋮ \\ v_{1}^{n} (x) & \dots & v_{n}^{n} (x) \end{matrix}]

Kdž najdeme partikulární řešení $y_{0} (x)$ rovnice $y^{'} = 𝐀 (x) y + b (x)$ , můžeme vyjádřit řešení alfové úlohy:

y (x; α) = Φ (x) α + y_{0} (x)

Teď budeme hledat správné $α$ :

𝐔 α + 𝐕 (Φ (b) α + y_{0} (b)) = c

(𝐔 + 𝐕 Φ (b)) α = c - 𝐕 y_{0} (b)

To je prostě lineární soustava, jejíž řešitelnost je daná povahou matice $𝐔 + 𝐕 Φ (b)$ .

Potíže metody střelby

Zjevná nevýhoda je v tom, že malé rozdíly na začátku mohou způsobit velké rozdíly na konci. Ukážeme si to na konkrétní úloze.

y^{″} - 100 y = 0, y (0) = 1, y (10) = 1

Metoda střelby na více cílů

To, co nám u metody střelby způsobilo potíže, byl příliš dlouhý interval. Tak co kdybychom si ho nějakým způsobem rozsekali?

Vezměme si jako typovou úlohu opět rovnici s nelineární vazbou:

y^{'} = f (x, y), r (y (a), y (b)) = 0, x \in (a, b)

Rozdělíme interval $⟨ a, b ⟩$ na $m$ částí, které nemusí být stejně dlouhé. Budeme opakovaně řešit počáteční úlohu na intervalech $⟨ x_{k - 1}, x_{k} ⟩, k \in \hat{m}$ s řešením $y (x; α_{k})$ . Úlohy budou vypadat takto:

y^{'} = f (x, y), y (x_{k - 1}) = α_{k}, x \in (x_{k - 1}, x_{k})

Zároveň je musíme nějak provázat mezi sebou. K tomu si přidáme podmínky:

y (x_{k}; α_{k + 1}) = α_{k}, k \in \hat{m - 1}

r (α_{1}, y (x_{m}; α_{m - 1})) = 0

Vzniklou soustavu budeme řešit – jak jinak – Newtonovou metodou. Máme funkci

F (α) = [\begin{matrix} y (x_{1}; α_{0}) - α_{1} \\ y (x_{2}; α_{1}) - α_{2} \\ ⋮ \\ y (x_{m - 1}; α_{m - 2}) - α_{m - 1} \\ r (α_{1}, y (x_{m}; α_{m - 1})) \end{matrix}]

Newtonova iterace bude mít tvar

α^{l + 1} = α^{l} - {(F^{'} (α^{l}))}^{- 1} F (α^{l})

Hodnotu $F (α^{l})$ určíme řešením jednotlivých úloh (třeba Runge-Kuttovou metodou). S derivací už to bude horší:

F^{'} (α) = {(\partial_{{(α_{k})}_{j}} F_{i} (α))}_{i \in \hat{(m n)}, k \in \hat{m}, j \in \hat{n}} = (\begin{matrix} 𝐆_{1} & - 𝐈 & \dots \\ 𝐆_{2} & - 𝐈 & \dots \\ 𝐆_{3} & ⋱ \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋱ & ⋮ \\ \dots & 𝐆_{m - 1} & - 𝐈 \\ 𝐀 & \dots & 𝐁 \end{matrix})

{(𝐆_{k})}_{i, j} ≔ \partial_{{(α_{k})}_{j}} y_{i} (x_{k}; α_{k}), 𝐀_{i, j} ≔ \partial_{ξ_{j}} r_{i} (α_{1}, y (x_{m}; α_{m})), 𝐁_{i, j} ≔ \sum_{s = 1}^{n} \partial_{η_{s}} r_{i} (α_{1}, y (x_{m}; α_{m})) \cdot \partial_{{(α_{m})}_{j}} y_{s} (x_{m}; α_{m})

Teď se pořádně nadechneme a odvodíme si úlohu ve variacích:

\partial_{x} y (x; α_{k}) = f (x, y (x; α_{k})), y (x_{k - 1}; α_{k}) = α_{k}, x \in (x_{k - 1}, x_{k})

\partial_{x} \partial_{{(α_{k})}_{j}} y_{i} (x; α_{k}) = \sum_{s = 1}^{n} \partial_{y_{j}} f_{i} (x, y (x; α_{k})) \cdot \partial_{{(α_{k})}_{j}} y_{j} (x; α_{k}), \partial_{{(α_{k})}_{j}} y_{i} (x_{k - 1}; α_{k}) = δ_{i, j}

Tyto úlohy můžeme bez záruky zkusit nahradit definicí derivace:

\partial_{{(α_{k})}_{j}} (x_{k}, α_{k}) \approx \frac{y_{i} (x_{k}; α_{k}) - y (x_{k}; α_{k} - ({(α_{k})}_{j} - {(α_{k - 1})}_{j}) e_{j})}{{(α_{k})}_{j} - {(α_{k - 1})}_{j}}

Metoda sítí

Alternativní názvy: metoda konečných diferencí, metoda konečných rozdílů

Metoda spočívá v tom, že aproximujeme derivace diferencemi v konečně mnoha bodech a dostaneme tím algebraicky řešitelnou soustavu rovnic.

Jako typovou úlohu si vezmeme naši oblíbenou samoadjungovanou rovnici:

- {(p (x) y^{'})}^{'} + q (x) y = f (x), y (a) = γ_{1}, y (b) = γ_{2}, x \in (a, b)

p, p^{'}, q, f \in 𝒞 (⟨ a, b ⟩), q (x) \geq 0, p (x) \geq c_{0} > 0, x \in ⟨ a, b ⟩

Nejprve diskretizujeme množinu $⟨ a, b ⟩$ tím, že na ni rozmístíme vnější a vnitřní síť uzlů (pro jednoduchost je rozmístíme rovnoměrně):

{\bar{ω}}_{h} ≔ {a + j h | j = 0, \dots, m}

ω_{h} ≔ {a + j h | j = 1, \dots, m - 1}

$m$ je počet dílů, $h$ je velikost oka sítě. Tyto dvě sítě se liší v tom, jestli obsahují hraniční uzly:

γ_{h} ≔ {\bar{ω}}_{h} - ω_{h} = {a, a + m h}

Diferenciální výrazy diskretizujeme pomocí Taylora (taky by to šlo pomocí Lagrange, ale Taylor je jednodušší na odvození a zapamatování). Nechť $\hat{f} \in 𝒞^{m} (⟨ 0, 1 ⟩)$ , potom

\hat{f} (1) = \sum_{k = 1}^{m - 1} \frac{{\hat{f}}^{k} (0)}{k!} + m \int_{0}^{1} \frac{s^{m - 1} {\hat{f}}^{m} (1 - s)}{m!} d s

(S tímto tvarem zbytku jsme se v MAN2 nesetkali, ale údajně se dá snadno odvodit.)

Pro naše účely si vezmeme $\hat{f} (t) ≔ y (x + t h), t \in ⟨ 0, 1 ⟩$ . Spočteme si

{\hat{f}}^{k} (t) = y^{k} (x + t h) h^{k}

Použitím vzorečku máme

y (x + h) = \sum_{k = 1}^{m - 1} \frac{y^{k} (x) h^{k}}{k!} + \underset{≕ R^{m} (x, h)}{\underset{⏟}{m \int_{0}^{1} \frac{s^{m - 1} y^{m} (x + (1 - s) h) h^{m}}{m!} d s}}

Všimněme si, že výraz $\frac{R^{m} (x, h)}{h^{m}}$ je omezený.

𝒪

ℴ

R^{m} (x, h) = 𝒪 (h^{m})

Diferenční náhrady nejběžnějších typů derivací

Náhrady $y^{'} (x)$

Pro $m = 2$ máme

y (x + h) = y (x) + y^{'} (x) h + R^{2} (x, h)

y^{'} (x) = \frac{y (x + h) - y (x)}{h} + 𝒪 (h)

Provedeme-li stejný postup s náhradou $h$ za $- h$ , dostaneme

y^{'} (x) = \frac{y (x) - y (x - h)}{h} + 𝒪 (h)

Pro $m = 3$ dostaneme

y (x \pm h) = y (x) \pm y^{'} (x) h + \frac{y^{″} (x) h^{2}}{2} + R^{3} (x, \pm h)

Odečteme-li tyto dvě rovnosti (pro $+$ a $-$ ) od sebe, dostaneme

y (x + h) - y (x - h) = 2 y^{'} (x) h + R^{3} (x, h) - R^{3} (x, - h)

y^{'} (x) = \frac{y (x + h) - y (x - h)}{2 h} + 𝒪 (h^{2})

Shrneme si všechny náhrady první derivace do tabulky:

Název	Náhrada	Chyba
dopředná diference	$y_{x} (x) ≔ \frac{y (x + h) - y (x)}{h}$	$𝒪 (h)$
zpětná diference	$y_{\bar{x}} (x) ≔ \frac{y (x) - y (x - h)}{h}$	$𝒪 (h)$
centrální diference	$y_{\overset{˚}{x}} (x) ≔ \frac{y (x + h) - y (x - h)}{2 h}$	$𝒪 (h^{2})$

Náhrada $y^{″} (x)$

Pro $m = 4$ máme

y (x \pm h) = y (x) \pm y^{'} (x) h + \frac{y^{″} (x) h^{2}}{2} \pm \frac{y^{‴} (x) h^{3}}{6} + R^{4} (x, \pm h)

Sečtením obou rovností dostáváme

y (x + h) + y (x - h) = 2 y (x) + y^{″} (x) h^{2} + R^{4} (x, h) + R^{4} (x, - h)

y^{″} (x) = \frac{y (x + h) - 2 y (x) + y (x - h)}{h^{2}} + 𝒪 (h^{2})

Zavedeme tedy centrální druhou diferenci:

y_{\bar{x} x} (x) ≔ \frac{y (x + h) - 2 y (x) + y (x - h)}{h^{2}}

Toto označení sice vypadá divně, ale dává smysl, protože můžeme snadno ověřit, že ${(y_{\bar{x}})}_{x} (x) = y_{\bar{x} x} (x)$ .

Náhrada $y^{4} (x)$

Podobným způsobem odvodíme centrální čtvrtou diferenci:

y^{4} (x) = y_{\bar{x} x \bar{x} x} (x) + 𝒪 (h^{2})

y_{\bar{x} x \bar{x} x} (x) ≔ \frac{y (x + 2 h) - 4 y (x + h) + 6 y (x) - 4 y (x - h) + y (x - 2 h)}{h^{4}}

V našem příkladu použijeme náhradu ${(p (x) y^{'})}^{'} \approx {(p y_{\bar{x}})}_{x} (x)$ (ale to není jediná možnost). Speciálně pro $p$ konstantní máme ${(p y_{\bar{x}})}_{x} (x) = p y_{\bar{x} x} (x)$ . Dokážeme si, že to dobře funguje i pro nekonstantní.

p, y

Nyní si můžeme sestavit algebraickou rovnici pro aproximaci řešení v uzlech sítě.

u : {\bar{ω}}_{h} \to ℝ

Diferenciální rovnici omezíme na uzly sítě ${\bar{ω}}_{h}$ :

- {(p y^{'})}^{'} (a + j h) + q (a + j h) y (a + j h) = f (a + j h), j = 1, \dots, m - 1

Provedeme-li náhradu konečnou diferencí a označíme si síťové funkce, dostaneme

- {({(p y_{\bar{x}})}_{x})}_{j} + 𝒪 (h^{α}) + q_{j} y_{j} = f_{j}

Se členem $𝒪 (h^{α})$ toho moc neuděláme, protože vůbec netušíme, co v něm je. Můžeme se spoléhat akorát na to, že když zvolíme dostatečně malé $h$ , tak bude taky dostatečně malý. Tudíž se na něj vykašleme. Ale tím pádem už nemáme přesné řešení, což zdůrazníme tím, že funkci $f$ přejmenujeme na $u$ . Rovnici poté můžeme vyjádřit ve vektorovém tvaru:

- {(p u_{\bar{x}})}_{x} + q u = f, u_{0} = γ_{1}, u_{m} = γ_{2}

Pokud si do toho dosadíme naše krásné vyjádření ${(p u_{\bar{x}})}_{x}$ , dostaneme soustavičku

\begin{align} u_{0} & = γ_{1} \\ - \frac{p_{1}}{h^{2}} u_{0} + \frac{p_{1} + p_{2}}{h^{2}} u_{1} - \frac{p_{2}}{h^{2}} u_{2} + q_{1} u_{1} & = f_{1} \\ - \frac{p_{2}}{h^{2}} u_{0} + \frac{p_{2} + p_{3}}{h^{2}} u_{2} - \frac{p_{3}}{h^{2}} u_{3} + q_{2} u_{2} & = f_{2} \\ ⋮ \\ - \frac{p_{j}}{h^{2}} u_{0} + \frac{p_{j} + p_{j + j}}{h^{2}} u_{j} - \frac{p_{j + j}}{h^{2}} u_{j + j} + q_{j} u_{j} & = f_{j} \\ ⋮ \\ - \frac{p_{m - 1}}{h^{2}} u_{0} + \frac{p_{m - 1} + p_{m}}{h^{2}} u_{m - 1} - \frac{p_{m}}{h^{2}} u_{m} + q_{m - 1} u_{m - 1} & = f_{m - 1} \\ u_{m} & = γ_{2} \end{align}

Vidíme, že vznikla soustava lineárních rovnic s tridiagonální maticí. Budeme ji řešit maticovou faktorizací. Kompletně si všechno přeznačíme, abychom měli dost písmenek:

u_{0} = ϰ_{1} u_{1} + μ_{1}

A_{j} u_{j - 1} - C_{j} u_{j} + B_{j} u_{j + 1} = - F_{j}, j = 1, \dots, m - 1

u_{m} = ϰ_{2} u_{m - 1} + μ_{2}

přičemž pro tuto konkrétní úlohu máme

ϰ_{1, 2} = 0, μ_{1, 2} = γ_{1, 2}

A_{j} = - \frac{p_{j}}{h^{2}}, C_{j} = - \frac{p_{j} + p_{j + 1}}{h^{2}} - q_{j}, B_{j} = - \frac{p_{j + 1}}{h^{2}}, F_{j} = - f_{j}, j = 1, \dots, m - 1

Teď to vypadá, že jame si tam $ϰ_{1, 2}$ zaváděli zbytečně, ale u jiné okrajové podmínky by se mohly objevit.

Řešení funguje tak, že zvolíme $α_{1} ≔ ϰ_{1}, β_{1} ≔ μ_{1}$ a budeme rekurentně počítat

α_{j + 1} ≔ \frac{B_{j}}{C_{j} - α_{j} A_{j}}, β_{j + 1} ≔ \frac{β_{j} A_{j} + F_{j}}{C_{j} - α_{j} A_{j}}, j = 1, \dots, m - 1

Potom si zpětně napočteme

u_{m} ≔ \frac{μ_{2} + ϰ_{2} β_{m}}{1 - α_{m} ϰ_{2}}

u_{j - 1} ≔ α_{j} u_{j} + β_{j}, j = m, \dots, 1

Tím jsme našli řešení. Teď vyvstává otázka, jak velké chyby jsme se dopustili. Chceme porovnat skutečné řešení $y$ s řešením v podobě síťové funkce $u$ , které jsme našli. Samozřejmě potíž je v tom, že $y$ je definovaná na $⟨ a, b ⟩$ , zatímco $u$ jen na ${\bar{ω}}_{h}$ . Zjevné řešení je prostě zúžit $y$ na ${\bar{ω}}_{h}$ . Druhá možnost je nějakým způsobem rozšířit $u$ na $⟨ a, b ⟩$ .

𝒫_{h} : 𝒞 (⟨ a, b ⟩) \to ℋ_{h}, 𝒫_{h} y ≔ y |_{{\bar{ω}}_{h}}

S_{h} : ℋ_{h} \to L^{2} ((a, b)), S_{h} u (x) ≔ u_{j}, x \in (a + j h - \frac{h}{2}, a + j h + \frac{h}{2})

Q_{h} : ℋ_{h} \to 𝒞 (⟨ a, b ⟩), Q_{h} u (x) ≔ u_{j - 1} + \frac{u_{j} - u_{j - 1}}{h} (x - a - (j - 1) h), x \in ⟨ a + (j - 1) h, a + j h ⟩

My budeme pro jednoduchost porovnávat jen pomocí operátoru restrikce. K tomu potřebujeme způsob, jak měřit síťové funkce.

h \in (0, 1 ⟩

Zapíšeme obecně přesné řešení:

L y = f, x \in (a, b), l y |_{x = a} = γ_{1}, l y |_{x = b} = γ_{2}

a řešení diferenční úlohy:

L_{h} u = f, x \in ω_{h}, l_{h} u |_{0} = γ_{1}, l_{h} u |_{m} = γ_{2}

Okrajové podmínky například pro $x = a$ můžou mít tvar:

Dirichletova: $l y |_{x = a} = y (a), l_{h} u |_{0} = u_{0}$
Neumannova: $l y |_{x = a} = p (a) y^{'} (a), l_{h} u |_{0} = p_{0} u_{x, 0}$
Newtonova: $l y |_{x = a} = α_{1} p (a) y^{'} (a) - β_{1} y (a), l_{h} u |_{0} = α_{1} p_{0} u_{x, 0} - β_{1} u_{0}$

Chyba aproximace diferenciálního operátoru

y \in 𝒞 (⟨ a, b ⟩)

L

ψ_{0} = α_{1} p_{o} {(𝒫_{h} y)}_{x, 0} - β_{1} {(𝒫_{h} y)}_{0} - 𝒫_{h} (α_{1} p (a) y^{'} (a) - β_{1} y (a)) = α_{1} p_{0} ({(𝒫_{h} y)}_{x, 0} - (𝒫_{h} y^{'} (a))) = 𝒪 (h)

{‖ \cdot ‖}_{h}

𝐀_{h} u = φ

h_{0} \in ℝ^{+}

u

u, v \in ℋ_{h}

u, v \in ℋ_{h}

p, u, v \in ℋ_{h}

p, u, v \in ℋ_{h}

{‖ u ‖}_{0, h} ≔ max_{j = 0, \dots, m} | u_{j} |

ℋ_{h}^{*} ≔ {u : {\bar{ω}}_{h} \to ℝ | u_{0} = u_{m} = 0}

ℋ_{h}^{*} ≔ {u : {\bar{ω}}_{h} \to ℝ | u_{0} = u_{m} = 0}

Poznámka První Sobolevova nerovnost platí jen v dimenzi 1, ale ta druhá platí v každé dimenzi (což si dokážeme později).

ℋ_{h}^{*} ≔ {u : {\bar{ω}}_{h} \to ℝ | u_{0} = u_{m} = 0}

Metoda energetických nerovností

Metoda energetických nerovností je způsob, jak určit stabilitu metody.

Vezměme naši oblíbenou úlohu:

- {(p y^{'})}^{'} + q y = f, x \in (a, b), y (a) = γ_{1}, y (b) = γ_{2}

K ní máme diferenční schéma:

- {(p u_{\bar{x}})}_{x} + q u = f, x \in ω_{h}, u_{0} = γ_{1}, u_{1} \in γ_{2}

Odečteme zprojektovanou diferenciální rovnici od diferenční rovnice:

- {(p u_{\bar{x}})}_{x} + q u + 𝒫_{h} ({(p y^{'})}^{'}) - 𝒫_{h} (q y) = 0, x \in ω_{h}, u_{0} = y (a), u_{m} = y (b)

Použijeme chybu diferenciálního operátoru:

ψ = L_{h} (𝒫_{h} y) - 𝒫_{h} (L y) = - {(p {(𝒫_{h} y)}_{\bar{x}})}_{x} + 𝒫_{h} ({(p y^{'})}^{'}), ψ_{0} = ψ_{m} = 0

- {(p u_{\bar{x}})}_{x} + {(p {(𝒫_{h} y)}_{\bar{x}})}_{x} + q (u - 𝒫_{h} y) + ψ = 0

- {(p {(u - 𝒫_{h} y)}_{\bar{x}})}_{x} + q (u - 𝒫_{h} y) = - ψ

Označíme-li $z ≔ u - 𝒫_{h} y$ , dostaneme

- {(p z_{\bar{x}})}_{x} + q z = - ψ, z_{0} = z_{m} = 0

Vynásobíme obě strany rovnice zprava $z$ v součinu ${(\cdot, \cdot)}_{h}$ :

- {({(p z_{\bar{x}})}_{x}, z)}_{h} + {(q z, z)}_{h} = - {(ψ, z)}_{h}

Použijeme bu bun seki-bun:

(p z_{\bar{x}}, z_{\bar{x}}] + {(q z, z)}_{h} = - {(ψ, z)}_{h}

Tomu se říká energetická rovnost, protože ve fyzice ty členy často mají rozměr energie.

U pravé strany aplikujeme Schwarzovu nerovnost $- {(ψ, z)}_{h} \leq {‖ ψ ‖}_{h} {‖ z ‖}_{h}$ . Na levé straně u druhého členu použijeme předpoklad $q \geq 0$ o původní diferenciální úloze:

{(q z, z)}_{h} = \sum_{j = 1}^{m - 1} h q_{j} z_{j}^{2} \geq 0

U prvního členu levé strany využijeme předpoklad $p > c_{0} \in ℝ^{+}$ o původní diferenciální úloze:

(p z_{\bar{x}}, z_{\bar{x}}] = \sum_{j = 1}^{m} h p_{j} {(z_{\bar{x}})}_{j}^{2} \geq c_{0} {‖ z_{\bar{x}} ⟧}^{2}

Dosazením těchto vyjádření do energetické rovnosti a použitím třetí Sobolevovy nerovnosti dostáváme

c_{0} {‖ z_{\bar{x}} ⟧}^{2} \leq {‖ ψ ‖}_{h} {‖ z ‖}_{h} \leq \frac{b - a}{\sqrt{8}} {‖ ψ ‖}_{h} ‖ z_{\bar{x}} ⟧

Vydělíme rovnici $c_{0} ‖ z_{\bar{x}} ⟧$ :

‖ z_{\bar{x}} ⟧ \leq \underset{konstanta stability}{\underset{⏟}{\frac{b - a}{c_{0} \sqrt{8}}}} {‖ ψ ‖}_{h}

Z první Sobolevovy nerovnosti potom plyne

{‖ z ‖}_{0, h} \leq \frac{\sqrt{b - a}}{2} ‖ z_{\bar{x}} ⟧ \leq \underset{M}{\underset{⏟}{\frac{\frac{{(b - a)}^{3}}{2}}{c_{0} \sqrt{32}}}} {‖ ψ ‖}_{h}

Použitím Laxovy věty dostáváme, že schéma je stabilní.

Metoda sítí pro další okrajové podmínky

Vezměme jednoduchou okrajovou úlohu s Newtonovou (alias Robinovou) podmínkou:

- y^{″} = f, x \in (a, b), - y^{'} (a) = β_{1} y (a) + γ_{1}, y^{'} (b) = β_{2} y (b) + γ_{2}

K ní si vyrobíme diferenční schéma:

- u_{\bar{x} x} = f, x \in ω_{h}, - {(u_{x})}_{0} = β_{1} u_{0} + γ_{1}, {(u_{\bar{x}})}_{m} = β_{2} u_{m} + γ_{2}

Vyrobíme si z toho soustavu lineárních rovnic:

\begin{align} - \frac{u_{1} - u_{0}}{h} - β_{1} u_{0} & = γ_{1} \\ - \frac{u_{2} - 2 u_{1} + u_{0}}{h^{2}} & = f_{1} \\ ⋮ \\ - \frac{u_{m} - 2 u_{m - 1} + u_{m - 2}}{h^{2}} & = f_{m - 1} \\ \frac{u_{m} - u_{m - 1}}{h} - β_{2} u_{m} = γ_{1} \end{align}

Ta má opět tridiagonální tvar, takže ji můžeme řešit metodou faktorizace. Ještě si z nejasného důvodu vynásobíme první a poslední rovnici $\frac{2}{h}$ a můžeme ji zapsat v maticovém tvaru $𝐀 u = φ$ :

𝐀 ≔ (\begin{matrix} - \frac{2 β_{1}}{h} + \frac{2}{h^{2}} & - \frac{2}{h^{2}} & \dots \\ - \frac{1}{h^{2}} & \frac{2}{h^{2}} & - \frac{1}{h^{2}} & \dots \\ ⋮ & ⋮ & ⋮ & ⋱ & - \frac{1}{h^{2}} & \frac{2}{h^{2}} & - \frac{1}{h^{2}} \\ ⋮ & ⋮ & ⋮ & ⋱ & - \frac{2}{h^{2}} & - \frac{2 β_{2}}{h} + \frac{2}{h^{2}} \end{matrix}), φ ≔ (\begin{matrix} \frac{2 γ_{1}}{h} \\ f_{1} \\ ⋮ \\ f_{m - 1} \\ \frac{2 γ_{2}}{h} \end{matrix})

Uvažujme množinu síťových funkcí $ℋ_{h}$ se skalárním součinem $[u, v] = \sum_{j = 1}^{m - 1} h u_{j} v_{j} - \frac{h}{2} (u_{0} v_{0} + u_{m} v_{m})$ . Jelikož nemáme čas, jenom si bez důkazu řekneme pár vět.

𝐀

C_{1} \in ℝ^{+}

𝒪 (\sqrt{h})

Metoda sítí pro eliptické parciální diferenciální rovnice

Nechť $Ω ≔ (0, L_{1}) \times (0, L_{2}), f : Ω \to ℝ, g : \partial Ω \to ℝ$ . Uvažujme typovou úlohu

- \partial_{1} (λ (x) \partial_{1} y) - \partial_{2} (λ (x) \partial_{2} y) + q (x) y = f (x), x \in Ω, y |_{\partial Ω} = g

s řešením $y : \bar{Ω} \to ℝ, y ≔ y (x)$ .

Budeme předpokládat $f \in 𝒞 (\bar{Ω}), g \in 𝒞 (\partial Ω), λ^{'}, q \in 𝒞 (\bar{Ω}), q \geq 0, λ \geq C_{0} \in ℝ^{+}$ . V Benešově světě samozřejmě může být funkce spojitá na množině, na které není definovaná.

Tentokrát už budeme potřebovat dvojrozměrnou síť:

{\bar{ω}}_{h} ≔ {[i h_{1}, j h_{2}] \in \bar{Ω} | i = 0, \dots, m_{1}, j \in 0, \dots, m_{2}}, h_{1} = \frac{L_{1}}{m_{1}}, h_{2} = \frac{L_{2}}{m_{2}}

Také si potřebujeme vytvořit dvourozměrné diferenční náhrady.

\partial_{i} (λ (x) \partial_{i} y) = {(λ y_{{\bar{x}}_{i}})}_{x_{i}} + 𝒪 (h_{i}^{α})

kde $α = 2$ pro konstantní $λ$ a $α = 1$ pro nekonstantní $λ$ .

Označíme-li $y_{i, j} ≔ y (i h_{1}, j h_{2})$ , potom

{({(λ y_{{\bar{x}}_{1}})}_{x_{1}})}_{i, j} = \frac{1}{h_{1}} (λ_{i + 1, j} \frac{y_{i + 1, j} - y_{i, j}}{h_{1}} - λ_{i, j} \frac{y_{i, j} - y_{i - 1, j}}{h_{1}})

{({(λ y_{{\bar{x}}_{2}})}_{x_{2}})}_{i, j} = \frac{1}{h_{2}} (λ_{i, j + 1} \frac{y_{i, j + 1} - y_{i, j}}{h_{2}} - λ_{i, j} \frac{y_{i, j} - y_{i, j - 1}}{h_{2}})

Teď už můžeme sepsat diferenční schéma:

- {(λ y_{{\bar{x}}_{1}})}_{x_{1}} - {(λ y_{{\bar{x}}_{2}})}_{x_{2}} + q u = f, x \in ω_{h}, u |_{γ_{h}} = g |_{γ_{h}}

Označíme-li ${\bar{\nabla}}_{h} u ≔ (u_{{\bar{x}}_{1}}, u_{{\bar{x}}_{2}})$ a $\nabla_{h} \cdot \vec{W} ≔ {(W_{x_{1}})}_{1} + {(W_{x_{2}})}_{2}$ , můžeme to zapsat v hezčím tvaru:

- \nabla_{h} \cdot (λ {\bar{\nabla}}_{h} u) + q u = f, x \in ω_{h}, u |_{γ_{h}} = g |_{γ_{h}}

Teď to zapíšeme po řádcích, čímž vznikne naprosto nádherná soustava rovnic:

\underset{E_{i, j}}{\underset{⏟}{(\frac{λ_{i + 1, j} + λ_{i, j}}{h_{1}^{2}} + \frac{λ_{i, j + 1} + λ_{i, j}}{h_{2}^{2}} + q_{i, j})}} u_{i, j} \underset{A_{i, j}}{\underset{⏟}{- \frac{λ_{i, j}}{h_{1}^{2}}}} u_{i - 1, j} \underset{C_{i, j}}{\underset{⏟}{- \frac{λ_{i + 1, j}}{h_{1}^{2}}}} u_{i + 1, j} \underset{B_{i, j}}{\underset{⏟}{- \frac{λ_{i, j}}{h_{1}^{2}}}} u_{i, j - 1} \underset{D_{i, j}}{\underset{⏟}{- \frac{λ_{i, j + 1}}{h_{1}^{2}}}} u_{i, j + 1} = \underset{F_{i, j}}{\underset{⏟}{f_{i, j}}}

Abychom tuto soustavu mohli řešit, musíme si ji přeindexovat, aby místo indexů $i, j$ byl jen jeden. To už tady nebudeme podrobněji rozebírat; některé počítačové řešiče to umí dělat automaticky. Na řešení je dobrá super-relaxační metoda, metoda konjugovaných gradientů nebo Jacobiho metoda. Také se dá vytvořit jakási analogie metody faktorizace, ale ta je prý dost pomalá.

Teď pojďme analyzovat stabilitu.

ℋ_{h}^{*} ≔ {v : {\bar{ω}}_{h} \to ℝ | v |_{γ_{h}} = 0}

ℋ_{h} ≔ {v : {\bar{ω}}_{h} \to ℝ | v |_{γ_{h}} = 0}

Pomocí těchto nerovností podobně jako u jednorozměrného případu dokážeme, že schéma je stabilní.

- 𝒫_{h} (\nabla \cdot (λ (x) \nabla y)) + \nabla_{h} \cdot (λ {\bar{\nabla}}_{h} u) + q (𝒫_{h} y - u) = 0, 𝒫_{h} y |_{γ_{h}} = u |_{γ_{h}}

Použijeme chybu aproximace:

ψ = 𝒫_{h} (L y) - L_{h} (𝒫_{h} y) = - 𝒫_{h} (\nabla \cdot (λ \nabla y)) + \nabla_{h} \cdot (λ {\bar{\nabla}}_{h} (𝒫_{h} y))

- \nabla_{h} \cdot (λ {\bar{\nabla}}_{h} (𝒫_{h} y)) + \nabla_{h} \cdot (λ {\bar{\nabla}}_{h} u) + q (𝒫_{h} y - u) = - ψ, (𝒫_{h} y - u) |_{γ_{h}} = 0

Označíme-li $z ≔ {\hat{P}}_{y} - u$ , úlohu přepíšeme do tvaru

- \nabla_{h} \cdot (λ {\bar{\nabla}}_{h} z) + q z = - ψ, z |_{γ_{h}} = 0

To vynásobíme $z$ , čímž získáme energetickou rovnost:

- {(\nabla_{h} \cdot (λ {\bar{\nabla}}_{h} z), z)}_{h} + {(q z, z)}_{h} = - {(ψ, z)}_{h}

Aplikujeme Greenovu formuli:

(λ {\bar{\nabla}}_{h} z, {\bar{\nabla}}_{h} z] + {(q z, z)}_{h} = - {(ψ, z)}_{h}

Využijeme předpokladu, že $\forall x \in \bar{Ω} : q (x) \geq 0, λ (x) \geq c_{0} > 0$ , čímž odhadneme

{(q z, z)}_{h} \geq 0

\begin{align} (λ {\bar{\nabla}}_{h} z, {\bar{\nabla}}_{h} z] & = (λ z_{{\bar{x}}_{1}}, z_{{\bar{x}}_{1}} ⌋ + (λ z_{{\bar{x}}_{2}}, z_{{\bar{x}}_{2}} ⌉ \\ = \sum_{j = 1}^{m_{2} - 1} \sum_{i = 1}^{m_{1}} h_{1} h_{2} λ_{i, j} {| {(z_{{\bar{x}}_{1}})}_{i, j} |}^{2} + \sum_{i = 1}^{m_{1} - 1} \sum_{j = 1}^{m_{2}} h_{1} h_{2} λ_{i, j} {| {(z_{{\bar{x}}_{2}})}_{i, j} |}^{2} \\ \geq c_{0} ({\bar{\nabla}}_{h} z, {\bar{\nabla}}_{h} z] = c_{0} {‖ {\bar{\nabla}}_{h} z ⟧}^{2} \end{align}

Z energetické rovnosti a druhé Sobolevovy nerovnosti plyne

c_{0} {‖ {\bar{\nabla}}_{h} z ⟧}^{2} \leq - {(ψ, z)}_{h} \leq {‖ ψ ‖}_{h} {‖ z ‖}_{h} \leq c_{2, 2} {‖ ψ ‖}_{h} ‖ {\bar{\nabla}}_{h} z ⟧

‖ {\bar{\nabla}}_{h} z ⟧ \leq \frac{c_{2, 2}}{c_{0}} {‖ ψ ‖}_{h}

Případně můžeme nějak podobně odvodit

{‖ z ‖}_{h} \leq c_{2, 2} ‖ {\bar{\nabla}}_{h} z ⟧ \leq \frac{c_{2, 2}^{2}}{c_{0}} {‖ ψ ‖}_{h}

Metoda sítí pro evoluční úlohy

Vezmeme si typovou úlohu:

\partial_{t} y = D \partial_{x, x} y + f (t, x), x \in (0, T) \times (a, b)

y |_{x = a} = γ_{1}, y |_{x = b} = γ_{2}

y |_{t = 0} = y_{ini} (x)

Diskretizujeme si časoprostor:

{\bar{ω}}_{h} ≔ {a + j h | j = 0, \dots, m, h ≔ \frac{b - a}{m}}

ω_{h} ≔ {a + j h | j = 1, \dots, m - 1, h ≔ \frac{b - a}{m}}

γ_{h} ≔ {\bar{ω}}_{h} - ω_{h}

t_{k} ≔ k \cdot τ, k = 0, \dots, N_{T}, τ ≔ \frac{T}{N_{T}}

Pro $ω \in 𝒞 (⟨ 0, T ⟩ \times ⟨ a, b ⟩)$ si zavedeme značení

ω_{j}^{k} ≔ ω (k τ, a + j h)

ω ≔ [ω_{0}^{k}, \dots, ω_{m}^{k}]

\hat{ω} ≔ [ω_{0}^{k + 1}, \dots, ω_{m}^{k + 1}]

Použijeme diferenční náhradu

\partial_{x, x} y |_{j}^{k} = {(y_{\bar{x} x})}_{j}^{k} + 𝒪 (h^{2})

\partial_{t} y |_{j}^{k} = {\begin{cases} {(y_{t})}_{j}^{k} + 𝒪 (τ) = \frac{y_{j}^{k + 1} - y_{j}^{k}}{τ} + 𝒪 (τ) \\ {(y_{\bar{t}})}_{j}^{k} + 𝒪 (τ) = \frac{y_{j}^{k} - y_{j}^{k - 1}}{τ} + 𝒪 (τ) \\ {(y_{\overset{˚}{t}})}_{j}^{k} + 𝒪 (τ^{2}) = \frac{y_{j}^{k + \frac{1}{2}} - y_{j}^{k - \frac{1}{2}}}{τ} + 𝒪 (τ^{2}) \end{cases}

Explicitní schéma

Použijeme dopřednou diferenci.

\frac{\hat{u} - u}{τ} = D u_{\bar{x} x} + f, {\hat{u}}_{0} = γ_{1}, {\hat{u}}_{m} = γ_{2}, u |_{j}^{0} = y_{ini} (a + j h)

Chyba aproximace je $𝒪 (h^{2} + τ)$ .

Jelikož můžeme přímo napočítat $\hat{u} = u + τ D u_{\bar{x} x} + τ f$ , toto schéma se snadno implementuje. Zapíšeme si ho bodově:

u_{j}^{k + 1} = u_{j}^{k} + τ D \frac{u_{j + 1}^{k} - 2 u_{j}^{k} + u_{j - 1}^{k}}{h^{2}} + τ f_{j}^{k}

u_{0}^{k + 1} = γ_{1}, u_{m}^{k + 1} = γ_{1}

u_{j}^{0} = y_{ini} (a + j h)

Označíme si matičku, abychom to mohli jednoduše zapsat:

𝐀 ≔ [\begin{matrix} 1 - \frac{2 τ D}{h^{2}} & \frac{τ d}{h^{2}} & 0 & \dots \\ \frac{τ d}{h^{2}} & ⋱ & ⋱ & ⋱ \\ 0 & ⋱ & ⋱ & ⋱ \end{matrix}]

u^{k + 1} = 𝐀 u^{k} + τ f^{k} = 𝐀 (𝐀 u^{k - 1} + τ f^{k - 1}) + τ f^{k} = 𝐀^{2} u^{k - 1} + τ 𝐀 f^{k - 1} + τ f^{k} = 𝐀^{3} u^{k - 2} + τ 𝐀^{2} f^{k - 2} + τ 𝐀 f^{k - 1} + τ f^{k}

Aby byla metoda stabilní, potřebujeme, aby vlastní čísla matice ležela v intervalu $(- 1, 1)$ . Z našich úvah o úloze vlastních čísel plyne, že $𝐀$ má takováto vlastní čísla:

{\bar{λ}}_{l} ≔ 1 - \frac{4 τ D}{h^{2}} sin {(\frac{π l}{2 m})}^{2}

Chceme tedy, aby platilo

- 1 < 10 \frac{4 τ D}{h^{2}} sin {(\frac{π l}{2 m})}^{2} < 1

\frac{4 τ D}{h^{2}} sin {(\frac{π l}{2 m})}^{2} < 2

\frac{4 τ D}{h^{2}} < 2

To je prý dost striktní požadavek, takže s takovouhle metodou to jde pomalu.

Implicitní schéma

Tentokrát použijeme zpětnou diferenci.

\frac{\hat{u} - u}{τ} = D {\hat{u}}_{\bar{x} x} + \hat{f}

Chyba aproximace je opět $𝒪 (h^{2} + τ)$ . Ale tentokrát je $\hat{u}$ vyjáďřeno jen implicitně, takže pro jeho zjistění budeme muset řešit soustavu rovnic:

\hat{u} - τ D {\hat{u}}_{\bar{x} x} = u + τ \hat{f}

Bodový zápis:

u_{j}^{k + 1} - τ D \frac{u_{j + 1}^{k + 1} - 2 u_{j}^{k + 1} + u_{j - 1}^{k + 1}}{h^{2}} = u_{j}^{k} + τ f_{j}^{k - 1}

u_{0}^{k + 1} = γ_{1}, u_{m}^{k + 1} = γ_{1}

u_{j}^{0} = y_{ini} (a + j h)

Přepis do maticového tvaru:

𝐁 ≔ [\begin{matrix} 1 + \frac{2 τ D}{h^{2}} & - \frac{τ D}{h^{2}} & 0 & \dots \\ - \frac{τ D}{h^{2}} & ⋱ & ⋱ & ⋱ \\ 0 & ⋱ & ⋱ & ⋱ \end{matrix}]

𝐁 u^{k + 1} = u^{k} + τ f^{k + 1}

U^{k + 1} = B^{- 1} (u^{k} + τ f^{k + 1}) = \dots = {(B^{- 1})}^{k + 1} u^{0} + \dots

Opět chceme, aby spektrum $𝐁$ bylo v intervalu $(- 1, 1)$ . Tenotokrát si ale bez důkazu řekneme, že tam bude vždy.

Crautovo-Nicolsonové schéma

V podstatě je to aritmetický průměr předchozích dvou schémat.

\frac{\hat{u} - u}{τ} = \frac{{\hat{u}}_{\bar{x} x} + u_{\bar{x} x}}{2} + \frac{\hat{f} + f}{2}

Chyba aproximace je $𝒪 (h^{2} + τ^{2})$ , což je lepší než u předchozích dvou (jelikož v Taylorovi se něco vzájemně požere).

Bodový zápis je extrémně krásný:

u_{j}^{k + 1} + τ \frac{D τ}{h^{2}} (u_{j + 1}^{k + 1} - 2 u_{j}^{k + 1} + u_{j - 1}^{k + 1}) = u_{j}^{k} + \frac{D τ}{2 h^{2}} (u_{j + 1}^{k} - 2 u_{j}^{k} + u_{j - 1}^{k}) + \frac{τ}{2} (f_{j}^{k + 1} + f_{j}^{k})

u_{0}^{k + 1} = γ_{1}, u_{m}^{k + 1} = γ_{1}

u_{j}^{0} = y_{ini} (a + j h)

Dokážeme pomocí metody energetických nerovností a spousty halucinogenů, že schéma je stabilní nepodmíněně.

\partial_{t} y (a + j h, (k + \frac{1}{2} τ)) - \frac{u_{j}^{k + 1} - u_{j}^{k}}{τ} = D \partial_{x, x} y_{j}^{k + \frac{1}{2}} - \frac{D}{2} ({(u_{\bar{x} x})}_{j}^{k + 1} + {(u_{\bar{x} x})}_{j}^{k}) + f_{j}^{k + \frac{1}{2}} - \frac{1}{2} (f_{j}^{k + 1} - f_{j}^{k})

ψ = 𝒫_{h}^{k + \frac{1}{2}} (\partial_{t} y - D \partial_{x, x} y) - \frac{y^{k + 1} - y^{k}}{τ} + \frac{D}{2} ({(y_{\bar{x} x})}^{k + 1} + {(y_{\bar{x} x})}^{k}) - f^{k + \frac{1}{2}} + \frac{1}{2} (f^{k + 1} + f^{k})

\frac{y_{j}^{k + 1} - y_{j}^{k}}{τ} - \frac{u_{j}^{k + 1} - u_{j}^{k}}{τ} - \frac{D}{2} ({(y_{\bar{x} x})}_{j}^{k + 1} + {(y_{\bar{x} x})}_{j}^{k}) + \frac{D}{2} ({(u_{\bar{x} x})}_{j}^{k + 1} + {(u_{\bar{x} x})}_{j}^{k}) = - ψ

y_{0}^{k + 1} = u_{0}^{k + 1}, y_{m}^{k + 1} = u_{m}^{k + 1}, y_{j}^{0} = u_{j}^{0}

Označíme-li jako obvykle $z ≔ 𝒫_{h} y - u$ , docela se nám to vyhezčí:

\underset{{(z_{\bar{t}})}_{j}^{k + 1}}{\underset{⏟}{\frac{z_{j}^{k + 1} - z_{j}^{k}}{τ}}} - \frac{D}{2} ({(z_{\bar{x} x})}_{j}^{k + 1} + {(z_{\bar{x} x})}_{j}^{k}) = - ψ_{j}^{k}, z_{0}^{k + 1} = 0, z_{m}^{k + 1} = 0, z_{j}^{0} = 0

Vynásobíme všechno zprava $τ {(z_{\bar{t}})}^{k + 1}$ :

τ {({(z_{\bar{t}})}^{k + 1}, {(z_{\bar{t}})}^{k + 1})}_{h} - \frac{D}{2} {({(z_{\bar{x} x})}^{k + 1} + {(z_{\bar{x} x})}^{k}, τ {(z_{\bar{t}})}^{k + 1})}_{h} = - {(ψ^{k}, τ {(z_{\bar{t}})}^{k + 1})}_{h}

Použijeme Greenovu větu:

τ {‖ {(z_{\bar{t}})}^{k + 1} ‖}_{h}^{2} + \frac{D}{2} \underset{{‖ {(z_{\bar{x}})}^{k + 1} ⟧}^{2} - {‖ {(z_{\bar{x}})}^{k} ⟧}^{2}}{\underset{⏟}{({(z_{\bar{x}})}^{k + 1} + {(z_{\bar{x}})}^{k}, {(z_{\bar{x}})}^{k + 1} - {(z_{\bar{x}})}^{k}]}} = - τ {(ψ^{k}, {(z_{\bar{t}})}^{k + 1})}_{h}

Podle Schwarzovy a Youngovy nerovnosti máme

- τ {(ψ^{k}, {(z_{\bar{t}})}^{k + 1})}_{h} \leq τ {‖ ψ^{k} ‖}_{h} {‖ {(z_{\bar{t}})}^{k + 1} ‖}_{h} \leq \frac{τ}{2} {‖ ψ^{k} ‖}_{h}^{2} + \frac{τ}{2} {‖ {(z_{\bar{t}})}^{k + 1} ‖}_{h}^{2}

Když to poskládáme dohromady, odečtením $\frac{τ}{2} {‖ {(z_{\bar{t}})}^{k + 1} ‖}_{h}^{2}$ od obou stran dostáváme

\frac{τ}{2} {‖ {(z_{\bar{t}})}^{k + 1} ‖}_{h}^{2} + \frac{D}{2} {‖ {(z_{\bar{x}})}^{k + 1} ⟧}^{2} - \frac{D}{2} {‖ {(z_{\bar{x}})}^{k} ⟧}^{2} \leq \frac{τ}{2} {‖ ψ^{k} ‖}_{h}^{2}

Z levé strany nerovnosti jistě můžeme bez újmy na platnosti smazat první člen, jelikož je kladný, a vše vynásobit dvěma, čímž dostaneme

D {‖ {(z_{\bar{x}})}^{k + 1} ⟧}^{2} \leq D {‖ {(z_{\bar{x}})}^{k} ⟧}^{2} + τ {‖ ψ^{k} ‖}_{h}^{2}

Iterováním této nerovnosti máme

D {‖ {(z_{\bar{x}})}^{k + 1} ⟧}^{2} \leq D {‖ {(z_{\bar{x}})}^{0} ⟧}^{2} + \sum_{j = 0}^{k} τ \underset{𝒪 (τ^{2} + h^{2})}{\underset{⏟}{{‖ ψ^{j} ‖}_{h}^{2}}} \overset{τ, h \to 0}{\to} 0

Z toho plyne, že metoda je stabilní.

Metoda přímek

To úplně nesouvisí s předchozími třemi schématy, ale je to jiný způsob, jak řešit tu samou úlohu.

Diskretizujeme úlohu v $(a, b)$ , tedy čas necháme spojitý. Tím vznikne systém obyčejných diferenciálních rovnic, který pak řešíme vhodným řešičem.

\partial_{x, x} y \approx y_{\bar{x} x} + 𝒪 (h^{2})

𝕦 = D u_{\bar{x} x} + f (t), u_{0} (t) = γ_{1}, u_{m} (t) = γ_{2}, u_{j} (0) = y_{ini} (a + j h)

Řešení této úlohy bude nějaká síťová funkce $u$ , která ale závisí na čase.

Bodově zapíšeme:

𝕦_{j} (t) = \frac{D}{h^{2}} (u_{j + 1} (t) - 2 u_{j} (t) + u_{j - 1} (t)) + f (t, a + j h), u_{0} (t) = γ_{1}, u_{m} (t) = γ_{2}, u_{j} (0) = y_{ini} (a + j h)

To je soustava $m - 1$ obyčejných diferenciálních rovnic, kterou můžeme řešit Eulerovou nebo Runge-Kuttovou metodou.

Numerická matematika 2

Obsah

Okrajové úlohy pro diferenciální rovnice

Typy běžných okrajových podmínek

Převod mezi tvary rovnic

Shrnutí poznatků o počátečních úlohách

Metoda střelby

Metoda střelby pro okrajovou podmínku nelineárního typu

Metoda střelby pro soustavu lineárních rovnic

Potíže metody střelby

Metoda střelby na více cílů

Metoda sítí

Diferenční náhrady nejběžnějších typů derivací

Náhrady $y^{'} (x)$

Náhrada $y^{″} (x)$

Náhrada $y^{4} (x)$

Chyba aproximace diferenciálního operátoru

Metoda energetických nerovností

Metoda sítí pro další okrajové podmínky

Metoda sítí pro eliptické parciální diferenciální rovnice

Metoda sítí pro evoluční úlohy

Explicitní schéma

Implicitní schéma

Crautovo-Nicolsonové schéma

Metoda přímek

Numerické řešení úloh s advekcí

Zákony zachování

Numerická matematika 2

Obsah

Okrajové úlohy pro diferenciální rovnice

Typy běžných okrajových podmínek

Převod mezi tvary rovnic

Shrnutí poznatků o počátečních úlohách

Metoda střelby

Metoda střelby pro okrajovou podmínku nelineárního typu

Metoda střelby pro soustavu lineárních rovnic

Potíže metody střelby

Metoda střelby na více cílů

Metoda sítí

Diferenční náhrady nejběžnějších typů derivací

Náhrady y′(x)

Náhrada y″(x)

Náhrada y4(x)

Chyba aproximace diferenciálního operátoru

Metoda energetických nerovností

Metoda sítí pro další okrajové podmínky

Metoda sítí pro eliptické parciální diferenciální rovnice

Metoda sítí pro evoluční úlohy

Explicitní schéma

Implicitní schéma

Crautovo-Nicolsonové schéma

Metoda přímek

Numerické řešení úloh s advekcí

Zákony zachování

Náhrady $y^{'} (x)$

Náhrada $y^{″} (x)$

Náhrada $y^{4} (x)$