Numerická matematika 2

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:

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 ${\overline{α}}_{1}, {\overline{α}}_{2}$ tak, aby $y (b; {\overline{α}}_{1}) < γ_{2}$ a $y (b; {\overline{α}}_{2}) > γ_{2}$ . Nastavíme $α_{1} \leftarrow {\overline{α}}_{1}, α_{2} \leftarrow {\overline{α}}_{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})

Φ (α) ≔ {‖ F (α) ‖}_{2}^{2} = \sum_{j = 1}^{n} {| F_{j} (α) |}^{2}

φ (λ) ≔ Φ (α_{k} + λ Δ α_{k})

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})}^{T} \cdot (F^{'} (α_{k} + λ Δ α_{k}) Δ α_{k})

φ^{'} (0) = 2 F {(α_{k})}^{T} \cdot (F^{'} (α_{k}) Δ α_{k}) = - 2 F {(α_{k})}^{T} \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 neorvnost 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 $ξ, η$ .

Příklad.

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

y_{2}^{'} = - 10 x y_{1} - 20 - x^{2}

y_{1} {(a)}^{2} - 2 = y_{2} {(b)}^{2} - 4 = 0

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

Dostaneme 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}

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{array}{ccc} v_{1}^{1} (x) & \dots & v_{n}^{1} (x) \\ ⋮ & ⋱ & ⋮ \\ v_{1}^{n} (x) & \dots & v_{n}^{n} (x) \end{array}]

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.

Příklad. Mějme diferenciální rovnici

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

Pro ilustraci ji nejprve vyřešíme analyticky. Máme fundamentální systém

y (x) = C_{1} e^{10 x} + C_{2} e^{- 10 x}

. Z okrajových podmínek určíme konstanty:

y (x) = \frac{1 - e^{- 100}}{e^{100} - e^{- 100}} e^{10 x} + \frac{e^{100} - 1}{e^{100} - e^{- 100}} e^{- 10 x}

Zkusme nyní použít metodu střelby. Budeme tedy řešit počáteční úlohu

y^{''} - 100 y = 0, y (0) = 1, y^{'} (0) = α

Tu ve skutečnosti taky vyřešíme analyticky, abychom získali obecnou představu, jak to bude dopadat pro různé hodnoty

α

. Řešení bude opět ve tvaru

y (x; α) = C_{1} e^{10 x} + C_{2} e^{- 10 x}

a dopočtením konstant dostaneme

y (x; α) = \frac{10 + α}{20} e^{10 x} + \frac{10 - α}{20} e^{- 10 x}

Dosazením druhé okrajové podmínky můžeme dopočíst správnou hodnotu

α

\frac{10 + α^{*}}{20} e^{100} + \frac{10 - α^{*}}{20} e^{- 100} = 1

α^{*} = \frac{20 - 10 e^{100} - 10 e^{- 100}}{e^{100} - e^{- 100}} = - 10 + 20 \frac{1 - e^{- 100}}{e^{100} - e^{- 100}}

Zlomek vpravo je téměř

0

, takže při numerickém výpočtu na počítači se toto číslo zaokrouhlí na

α_{0} = - 10

. Podívejme se, kolik nám v takovém případě vyjde druhá okrajová podmínka:

y (10; α_{0}) = \frac{10 - 10}{20} e^{100} + \frac{10 + 10}{20} e^{- 100} = e^{- 100}

To je téměř

0

, takže k

1

to má opravdu hodně daleko. Tudíž prakticky pro tento příklad nemá metoda střelby šanci se trefit.

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{array}{c} 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{array}]

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{array}{cccccc} 𝐆_{1} & - 𝐈 & \dots \\ 𝐆_{2} & - 𝐈 & \dots \\ 𝐆_{3} & ⋱ \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋱ & ⋮ \\ \dots & 𝐆_{m - 1} & - 𝐈 \\ 𝐀 & \dots & 𝐁 \end{array})

{(𝐆_{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ě):

{\overline{ω}}_{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} ≔ {\overline{ω}}_{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ý.

Definice (Landaův symbol $𝒪$ ). Pokud funkce

g

na okolí nuly

H_{0}

splňuje

(\exists K \in ℝ^{+}) (\exists α \in ℝ_{0}^{+}) (\forall x \in H_{0}) (| \frac{g (x)}{x^{α}} | \leq K)

potom říkáme, že

g

se na

H_{0}

chová jako

𝒪 (x^{α})

a píšeme

g (x) = 𝒪 (x^{α})

Definice (Landaův symbol $ℴ$ ). Pokud funkce

g

na okolí nuly

H_{0}

pro nějaké

α \in ℝ_{0}^{+}

splňuje

\lim_{x \to 0} \frac{g (x)}{x^{α}} = 0

potom říkáme, že

g

se na

H_{0}

chová jako

ℴ (x^{α})

a píšeme

g (x) = ℴ (x^{α})

Poznámka.

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_{\overline{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_{\overline{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_{\overline{x}})}_{x} (x) = y_{\overline{x} x} (x)$ .

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

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

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

y_{\overline{x} x \overline{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_{\overline{x}})}_{x} (x)$ (ale to není jediná možnost). Speciálně pro $p$ konstantní máme ${(p y_{\overline{x}})}_{x} (x) = p y_{\overline{x} x} (x)$ . Dokážeme si, že to dobře funguje i pro nekonstantní.

Věta. Nechť

p, y

jsou dostatečně regulární, konkrétně

p \in 𝒞^{3}, y \in 𝒞^{4}

, potom

{(p y^{'})}^{'} (x) = {(p y_{\overline{x}})}_{x} (x) + 𝒪 (h)

Důkaz. Použitím Taylora dostáváme

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

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

Jelikož výraz

{\tilde{R}}^{(3)} (x, h)

obsahuje třetí derivaci

p

a čtvrtou derivaci

y

, k jeho asymptotickému omezení požadujeme

p \in 𝒞^{3}, y \in 𝒞^{4}

. Vyjádříme si

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

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

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

Dosadíme, čímž vznikne extrémně krásný výraz:

{(p y^{'})}^{'} (x) = \frac{1}{h} (p (x + h) \frac{y (x + h) - y (x)}{h} - p (x) \frac{y (x) - y (x - h)}{h}) - \underset{𝒪 (h)}{\underset{⏟}{\frac{{(p y^{'})}^{''} (x) h}{2}}} - \underset{𝒪 (h^{2})}{\underset{⏟}{\frac{{\tilde{R}}^{(3)} (x, h)}{h}}} + \frac{1}{h} (p (x + h) (\underset{𝒪 (h)}{\underset{⏟}{\frac{y^{''} (x + h) h}{2}}} + \underset{𝒪 (h^{2})}{\underset{⏟}{\frac{R^{(3)} (x + h, - h)}{h}}}) - p (x) (\underset{𝒪 (h)}{\underset{⏟}{\frac{y^{''} (x) h}{2}}} + \underset{𝒪 (h^{2})}{\underset{⏟}{\frac{R^{(3)} (x, - h)}{h}}}))

Teď máme trochu problém, protože výrazy třídy

𝒪 (h)

v poslední hnusné závorce po vydělení

h

dávají

𝒪 (1)

. Budeme si muset poradit tím, že jejich rozdíl umlátíme pomocí věty o přírůstku funkce.

\frac{1}{h} (\dots - \dots) = \frac{p (x + h) y^{''} (x + h) - p (x) y^{''} (x)}{2} = \frac{{(p y^{''})}^{'} (ξ) h}{2} = 𝒪 (h), ξ \in (x, x + h)

Tím už jsme všechny nechtěné výrazy asymptoticky omezili jako

𝒪 (h)

, čímž dostáváme znění věty.

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

Definice. Síťová funkce je zobrazení

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

. Značíme

u (a + j h) ≕ u_{j}

. Potom

(u_{0}, \dots, u_{m})

je vektor síťové funkce.

Diferenciální rovnici omezíme na uzly sítě ${\overline{ω}}_{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_{\overline{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_{\overline{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_{\overline{x}})}_{x}$ , dostaneme soustavičku

\begin{aligned} 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{aligned}

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 ${\overline{ω}}_{h}$ . Zjevné řešení je prostě zúžit $y$ na ${\overline{ω}}_{h}$ . Druhá možnost je nějakým způsobem rozšířit $u$ na $⟨ a, b ⟩$ .

Definice. Nechť

ℋ_{h}

je množina všech síťových funkcí na

{\overline{ω}}_{h}

. Potom operátor restrikce je zobrazení

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

Definice. Nechť

ℋ_{h}

je množina všech síťových funkcí na

{\overline{ω}}_{h}

. Potom operátor po částech konstantního rozšíření je zobrazení

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

(intuitivně: kolem každého bodu uděláme rovnou plošku)

Definice. Nechť

ℋ_{h}

je množina všech síťových funkcí na

{\overline{ω}}_{h}

. Potom operátor po částech lineárního rozšíření je zobrazení

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 ⟩

(intuitivně: body pospojujeme do lomené čáry)

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

Definice. Nechť pro každé

h \in (0,1 ⟩

{‖ \cdot ‖}_{h}

norma na

ℋ_{h}

. Tyto normy jsou konzistentní, pokud existuje norma

{‖ \cdot ‖}_{0}

𝒞 (⟨ a, b ⟩)

taková, že pro všechna

y \in C (⟨ a, b ⟩)

\lim_{h \to 0^{+}} {‖ 𝒫_{h} y ‖}_{h} = {‖ y ‖}_{0}

Příklad. Pro

u \in ℋ_{h}

označme

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

a vezměme obyčejnou maximovou normu:

{‖ y ‖}_{0} ≔ \max {| y (x) | | x \in ⟨ a, b ⟩}

Potom

{‖ \cdot ‖}_{0, h}

jsou konzistentní normy konvergující k

{‖ \cdot ‖}_{0}

Příklad. Pro

u \in ℋ_{h}

můžeme zkusit vzít eukleidovskou normu

{‖ u ‖}_{2^{'}, h} ≔ \sqrt{\sum_{j = 0}^{m} u_{j}^{2}}

Pokud ale do této normy budeme strkat například konstantní funkce s čím dál menším

h

, půjde to k nekonečnu. Tudíž tato norma není konzistentní. Spravíme to tím, že si ji upravíme:

{‖ u ‖}_{2^{'}, h} ≔ \sqrt{\sum_{j = 1}^{m - 1} h u_{j}^{2} + \frac{h}{2} (u_{0}^{2} + u_{m}^{2})}

Když si tohle nakreslíme, jedná se vlastně o Riemannův integrální součet, tudíž to bude konvergovat k normě

{‖ \cdot ‖}_{2}

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

Definice. Nechť

y \in 𝒞 (⟨ a, b ⟩)

. Chyba aproximace diferenciálního operátoru

L

uvnitř

(a, b)

l

na okrajích je síťová funkce

ψ : {\overline{ω}}_{h} \to ℝ

daná výrazy

ψ ≔ L_{h} (𝒫_{h} y) - 𝒫_{h} (L y), x \in ω_{h}

ψ_{0} ≔ l_{h} (𝒫_{h} y) |_{0} - (L y) |_{0}

ψ_{m} ≔ l_{h} (𝒫_{h} y) |_{m} - (L y) |_{m}

Poznámka. Pro konkrétní tvary

L

l

výše máme v uzlech

ω_{h}

ψ = - {(p {(𝒫_{h} y)}_{\overline{x}})}_{x} + q 𝒫_{h} y - 𝒫_{h} (- {(p y^{'})}^{'} + q y) = - {(p {(𝒫_{h} y)}_{\overline{x}})}_{x} + 𝒫_{h} ({(p y^{'})}^{'}) = 𝒪 (h^{α})

kde

α = 2

pro

p

konstantní a

α = 1

pro

p

nekonstantní.

Příklad. Pro Newtonovu okrajovou podmínku:

ψ_{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)

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

Definice. Nechť

‖ \cdot ‖

je konzistentní norma.

u

konverguje k

y

, pokud

\lim_{h \to 0^{+}} {‖ u - 𝒫_{h} y ‖}_{h} = 0

Pokud existuje

α \in ℝ^{+}

takové, že

‖ u - 𝒫_{h} y ‖ = 𝒪 (h^{α})

, potom

u

konverguje s přesností řádu

α

y

Poznámka. Diferenční úlohu lze zapsat jako

𝐀 u = φ

, kde

𝐀_{h}

je matice závisející na

{\overline{ω}}_{h}

. Řešme tu samou s jinou pravou stranou:

𝐀 v = ψ

. Obě řešení můžeme porovnat jako

u - v = 𝐀_{h}^{- 1} (φ - ψ)

Aby byla úloha stabilní, požadujeme, aby malá změna

φ \to ψ

vedla k malé změně

u \to v

bez ohledu na výběr sítě. Chceme tedy najít konstantu

M \in ℝ^{+}

nezávislou na

{\overline{ω}}_{h}

a vhodné konzistentní normy

{‖ \cdot ‖}_{1, h}, {‖ \cdot ‖}_{2, h}

takové, že

{‖ u - v ‖}_{1, h} \leq M {‖ φ - ψ ‖}_{2, h}

Definice. Diferenční schéma je korektní, pokud existuje

h_{0} \in ℝ^{+}

splňující

pro každé $h \in (0, h_{0}]$ existuje množina $Φ \subset ℝ^{m + 1}$ taková, že pro všechna $φ \in Φ$ má úloha $𝐀_{h} u = φ$ jediné řešení (neboli $𝐀_{h}$ je regulární)
existuje $M \in ℝ^{+}$ nezávislé na $h$ a konzistentní normy ${‖ \cdot ‖}_{1, h}, {‖ \cdot ‖}_{2, h}$ takové, že $\forall h \in (0, h_{0}), φ, ψ \in Φ : {‖ 𝐀_{h}^{- 1} φ - 𝐀_{h}^{- 1} ψ ‖}_{1, h} \leq M ‖ φ - ψ ‖$

Toto

M

se nazývá konstanta stability.

Věta (Laxova). Nechť diferenční schéma je korektní a jeho chyba aproximace v nějaké konzistentní normě konverguje k nule. Potom řešení

u

diferenční úlohy konverguje k přesnému řešení

y

diferenciální rovnice.

Důkaz. Odečteme diferenční úlohu od diferenciální úlohy zúžené na síť:

𝒫_{h} (L y) - L_{h} u = 0, 𝒫_{h} (l y) |_{0} - l_{h} u |_{0} = 0, 𝒫_{h} (l y) |_{m} - l_{h} u |_{m} = 0

Použijeme definici chyby aproximace:

L_{h} (𝒫_{h} y) - ψ - L_{h} u = 0, l_{h} (𝒫_{h} y) |_{0} - ψ_{0} - l_{h} u |_{0} = 0, l_{h} (𝒫_{h} y) |_{m} - ψ_{m} - l_{h} u |_{m} = 0

Z linearity

L_{h}, l_{h}

to můžeme přepsat do tvaru

L_{h} (𝒫_{h} y - u) = ψ, l_{h} (𝒫_{h} y - u) |_{0} = ψ_{0}, l_{h} (𝒫_{h} y - u) |_{0} = ψ_{0}

Označíme

z ≔ 𝒫_{h} y - u

L_{h} z = ψ, l_{h} z |_{0} = ψ_{0}, l_{h} z |_{0} = ψ_{0}

Z předpokladu korektnosti máme

{‖ z ‖}_{1, h} \leq M {‖ ψ ‖}_{2, h}

Z předpokladu konvergence je

\lim_{h \to 0^{+}} {‖ ψ ‖}_{2, h} = 0

, takže i

\lim_{h \to 0^{+}} {‖ 𝒫_{h} y - u ‖}_{1, h} = 0

Definice (skalární součin síťových funkcí). Nechť

ℋ_{h}

je množina všech síťových funkcí na

{\overline{ω}}_{h}

. Pro

u, v \in ℋ_{h}

si zavedeme extrémně prokleté značení

{(u, v)}_{h} ≔ \sum_{j = 1}^{m - 1} h u_{j} v_{j}, {‖ u ‖}_{h} ≔ \sqrt{{(u, u)}_{h}}

(u, v] ≔ \sum_{j = 1}^{m} h u_{j} v_{j}, ‖ u ⟧ ≔ \sqrt{(u, u]}

[u, v) ≔ \sum_{j = 0}^{m - 1} h u_{j} v_{j}, ⟦ u ‖ ≔ \sqrt{[u, u)}

[u, v] ≔ {(u, v)}_{h} + \frac{h}{2} (u_{0} v_{0} + u_{m} v_{m}), ⟦ u ⟧ ≔ \sqrt{[u, u]}

To, že u toho prvního je index

_{h}

a u ostatních ne, naneštěstí není překlep.

Věta (bu bun seki-bun). Pro

u, v \in ℋ_{h}

platí

{(u_{x}, v)}_{h} = u_{m} v_{m} - u_{1} v_{0} - (u, v_{\overline{x}}]

{(u_{\overline{x}}, v)}_{h} = u_{m - 1} v_{m} - u_{0} v_{0} - [u, v_{x})

Důkaz.

\begin{aligned} {(u_{x}, v)}_{h} & = \sum_{j = 1}^{m - 1} h (\frac{u_{j + 1} - u_{j}}{h}) v_{j} \\ = \sum_{j = 1}^{m - 1} u_{j + 1} v_{j} - \sum_{j = 1}^{m - 1} u_{j} v_{j} \\ = \sum_{j = 2}^{m} u_{j} v_{j - 1} - \sum_{j = 1}^{m - 1} u_{j} v_{j} \\ = \sum_{j = 1}^{m} u_{j} (v_{j - 1} - v_{j}) - u_{1} v_{0} + u_{m} v_{m} \\ = - \sum_{j = 0}^{m - 1} h u_{j} \frac{v_{j + 1} - v_{j}}{h} - u_{1} v_{0} + u_{m} v_{m} \\ = - (u, v_{\overline{x}}] - u_{1} v_{0} + u_{m} v_{m} \end{aligned}

\begin{aligned} {(u_{\overline{x}}, v)}_{h} & = \sum_{j = 1}^{m - 1} h (\frac{u_{j} - u_{j - 1}}{h}) v_{j} \\ = \sum_{j = 1}^{m - 1} u_{j} v_{j} - \sum_{j = 1}^{m - 1} u_{j - 1} v_{j} \\ = \sum_{j = 1}^{m - 1} u_{j} v_{j} - \sum_{j = 0}^{m - 2} u_{j} v_{j + 1} \\ = \sum_{j = 0}^{m - 1} u_{j} (v_{j} - v_{j + 1}) - u_{0} v_{0} + u_{m - 1} v_{m} \\ = - \sum_{j = 0}^{m - 1} h u_{j} \frac{v_{j + 1} - v_{j}}{h} - u_{0} v_{0} + u_{m - 1} v_{m} \\ = - [u, v_{x}) - u_{0} v_{0} + u_{m - 1} v_{m} \end{aligned}

Věta (první Greenova formule). Nechť

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

. Potom

{({(p u_{\overline{x}})}_{x}, v)}_{h} = - (p u_{\overline{x}}, v_{\overline{x}}] + p_{m} {(u_{\overline{x}})}_{m} v_{m} - p_{1} {(u_{\overline{x}})}_{1} v_{0}

Důkaz. Použijeme bu bun seki-bun s

u ≔ p u_{\overline{x}}, v ≔ v

Věta (druhá Greenova formule). Nechť

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

. Potom

{({(p u_{\overline{x}})}_{x}, v)}_{h} - {(u, {(p v_{\overline{x}})}_{x})}_{h} = p_{m} ({(u_{\overline{x}})}_{m} v_{m} - u_{m} {(v_{\overline{x}})}_{m}) - p_{1} ({(u_{\overline{x}})}_{1} v_{0} - u_{0} {(v_{\overline{x}})}_{1})

Důkaz. Použijeme první Greenovu formuli, nejdříve normálně a poté s prohozením

u \leftrightarrow v

. Potom jen upravíme to, co vyjde.

Poznámka. Budeme používat normy

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

{‖ u ‖}_{h}, ‖ u ⟧

definované výše.

Lemma (první Sobolevova nerovnost). Definujme

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

Potom existuje

C_{1} (a, b) \in ℝ^{+}

takové, že pro všechna

u \in ℋ_{h}

platí

{‖ u ‖}_{0, h} \leq C_{1} (a, b) ‖ u_{\overline{x}} ⟧

Důkaz. Nejprve pro ilustraci dokažme podobné tvrzení pro

f \in 𝒞 (⟨ a, b ⟩), f (a) = f (b) = 0

. Pro každé

x \in ⟨ a, b ⟩

platí

f (x) - f (a) = \int_{a}^{b} f^{'}

Podle Schwarzovy nerovnosti máme

| f (x) | \leq \int_{a}^{b} | f^{'} | \leq \sqrt{\int_{a}^{b} 1^{2}} \sqrt{\int_{a}^{b} {| f^{'} |}^{2}} \leq \sqrt{b - a} {‖ f^{'} ‖}_{2}

Nyní vezměme

u \in ℋ_{h}

. Potom

\sum_{j = 1}^{k} h {(u_{\overline{x}})}_{j} = \sum_{j = 1}^{k} h \frac{u_{j} - u_{j - 1}}{h} = \sum_{j = 1}^{k} u_{j} - \sum_{j = 1}^{k} u_{j - 1} = u_{k} - u_{0}

a podobně

\sum_{j = k + 1}^{m} h {(u_{\overline{x}})}_{j} = \sum_{j = k + 1}^{m} h \frac{u_{j} - u_{j - 1}}{h} = u_{m} - u_{k}

Jelikož jsme předpokládali

u_{0} = u_{m} = 0

, můžeme si z těchto vztahů vyjádřit

u_{k}

dvěma způsoby:

u_{k} = \sum_{j = 1}^{k} h {(u_{\overline{x}})}_{j} = - \sum_{j = k + 1}^{m} h {(u_{\overline{x}})}_{j}

Z toho zjevně pro libovolné

λ \in [0,1]

platí

u_{k}^{2} = λ {(\sum_{j = 1}^{k} h {(u_{\overline{x}})}_{j})}^{2} + (1 - λ) {(- \sum_{j = k + 1}^{m} h {(u_{\overline{x}})}_{j})}^{2}

Aplikujeme Schwarzovu nerovnost:

u_{k}^{2} \leq λ \underset{k h}{\underset{⏟}{(\sum_{j = 1}^{k} h)}} (\sum_{j = 1}^{k} h {(u_{\overline{x}})}_{j}^{2}) + (1 - λ) \underset{(m - k) h}{\underset{⏟}{(- \sum_{j = k + 1}^{m} h)}} (- \sum_{j = k + 1}^{m} h {(u_{\overline{x}})}_{j}^{2})

Zvolíme-li speciálně

λ ≔ 1 - \frac{k}{m}

, potom

λ k h = (1 - λ) (m - k) h

, takže

u_{k}^{2} \leq (1 - \frac{k}{m}) k h \underset{{‖ u_{\overline{x}} ⟧}^{2}}{\underset{⏟}{\sum_{j = 1}^{m} h {(u_{\overline{x}})}_{j}^{2}}} = \underset{≕ g (\frac{k}{m})}{\underset{⏟}{(1 - \frac{k}{m}) \frac{k}{m}}} (b - a) {‖ u_{\overline{x}} ⟧}^{2}

Snadno zjistíme, že funkce

g (x) = (1 - x) x

má maximum

\frac{1}{4}

, takže můžeme odhadnout

u_{k} \leq \frac{b - a}{4} {‖ u_{\overline{x}} ⟧}^{2}

. Z toho plyne tvrzení věty, kde volíme

C_{1} (a, b) ≔ \frac{b - a}{4}

Lemma (druhá Sobolevova nerovnost). Definujme

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

Potom existuje

C_{2} (a, b) \in ℝ^{+}

takové, že pro všechna

u \in ℋ_{h}

platí

{‖ u ‖}_{h} \leq C_{2} (a, b) ‖ u_{\overline{x}} ⟧

Důkaz. Použitím první Sobolevovy nerovnosti:

{‖ u ‖}_{h}^{2} = \sum_{j = 1}^{m - 1} h {| u_{j} |}^{2} \leq \sum_{j = 1}^{m - 1} h {‖ u ‖}_{0, h}^{2} \leq \sum_{j = 1}^{m - 1} h \frac{b - a}{4} {‖ u_{\overline{x}} ⟧}^{2} = (m - 1) h \frac{b - a}{4} ‖ u_{\overline{x}} ⟧ < \frac{{(b - a)}^{2}}{4} ‖ u_{\overline{x}} ⟧

Stačí zvolit

C_{2} (a, b) ≔ \frac{b - a}{2}

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

Lemma (třetí Sobolevova nerovnost). Definujme

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

Potom pro všechna

u \in ℋ_{h}

platí

\frac{h^{2}}{4} {‖ u_{\overline{x}} ⟧}^{2} \leq {‖ u ‖}_{h}^{2} \leq \frac{(b - a)}{8} {‖ u_{\overline{x}} ⟧}^{2}

Důkaz. Vyřešme si jen tak pro inspiraci diferenciální úlohu na vlastní čísla:

- y^{''} - λ y = 0, x \in (a, b), y (a) = y (b) = 0

Na základě moudré rady od mistra Beneše budeme hledat řešení ve tvaru

y (x) = \sin (α \cdot (x - a)), α \in ℝ

. Ke splnění okrajových podmínek musí být

α = α_{k} ≔ \frac{k π}{b - a}, k \in ℤ

Dosazením dopočteme

y (x) = y_{k} (x) ≔ \sin (\frac{k π}{b - a} (x - a)), λ = λ_{k} ≔ {(\frac{k π}{b - a})}^{2}

Takže úloha jde řešit pro

λ

v tomto tvaru. Přitom s

k = 0

je rovnice nudná a pro

- k

to vyjde stejně jako pro

k

, takže stačí brát

k \in ℕ

. Teď pojďme řešit podobnou diferenční úlohu:

- u_{\overline{x} x} - λ u = 0, x \in ω_{h}, u_{0} = u_{m} = 0

Řešení hledáme ve tvaru

u_{j} = \sin α j

. Okrajové podmínky vynutí

α = α_{k} ≔ \frac{k π}{m}, k \in ℤ

. Definujeme si tedy

{(u_{k})}_{j} ≔ \sin \frac{k π}{m} j

a s použitím vzorečku

\sin (A + B) - \sin (A - B) = 2 \cos A \sin B

a nějakého dalšího dostaneme

\begin{aligned} {({(u_{k})}_{\overline{x} x})}_{j} & = \frac{{(u_{k})}_{j + 1} - 2 {(u_{k})}_{j} + {(u_{k})}_{j - 1}}{h^{2}} \\ = \frac{\sin \frac{k π}{m} (j + 1) - 2 \sin \frac{k π}{m} j + \sin \frac{k π}{m} (j - 1)}{h^{2}} \\ = \frac{(\sin \frac{k π}{m} (j + 1) - \sin \frac{k π}{m} j) - (\sin \frac{k π}{m} j - \sin \frac{k π}{m} (j - 1))}{h^{2}} \\ = \frac{2 \cos \frac{k π}{m} (j + \frac{1}{2}) \sin \frac{k π}{2 m} - 2 \cos \frac{k π}{m} (j - \frac{1}{2}) \sin \frac{k π}{2 m}}{h^{2}} \\ = \frac{2 \sin \frac{k π}{2 m} (- 2 \sin \frac{k π}{m} j \sin \frac{k π}{2 m})}{h^{2}} \\ = - \underset{≕ λ_{k}}{\underset{⏟}{\frac{4}{h^{2}} {(\sin \frac{k π}{2 m})}^{2}}} \underset{{(u_{k})}_{j}}{\underset{⏟}{\sin \frac{k π}{m} j}} \end{aligned}

Zjevně opět stačí

k \in ℕ

. Ale pořád je divné, že máme nekonečně mnoho řešení, když řešíme konečnou úlohu. (Nebo aspoň Benešovi je to divné; mně by to možná bylo divné, kdybych měl tušení, proč tohle děláme.) Dosadíme si do řešení

k ≔ m

k ≔ m + 1

{(u_{m})}_{j} = \sin \frac{m π}{m} j = \sin π j = 0

{(u_{m + 1})}_{j} = \sin \frac{(m + 1) π}{m} j = \sin (π j + \frac{π}{m} j) = - \sin (π j - \frac{π}{m} j) = \sin \frac{(m - 1) π}{m} j = - {(u_{m - 1})}_{j}

λ_{m + 1} = \frac{4}{h^{2}} {(\sin \frac{(m + 1) π}{2 m})}^{2} = \frac{4}{h^{2}} {(\sin (\frac{π}{2} + \frac{π}{2 m}))}^{2} = \frac{4}{h^{2}} {(\sin (\frac{π}{2} - \frac{π}{2 m}))}^{2} = \frac{4}{h^{2}} {(\sin \frac{(m - 1) π}{2 m})}^{2} = λ_{m - 1}

Analogicky pro každé

l = 1, \dots, m - 1

u_{m + l} = - u_{m - l}, λ_{m + l} = λ_{m - l}

. Takže jediná řešení jsou ve tvaru

{(u_{k})}_{j} = \sin \frac{k π}{m} j, λ_{k} ≔ \frac{4}{h^{2}} {(\sin \frac{k π}{2 m})}^{2}

Dokážeme, že síťové funkce

u_{1}, \dots, u_{m - 1}

tvoří ortogonální (se součinem

{(\cdot, \cdot)}_{h}

) bázi

ℋ_{h}

. K tomu nám stačí dokázat, že jsou ortogonální, protože je jich správný počet. Dosadíme nějaké

u_{k}

do rovnice a vynásobíme ji nějakým

u_{l}

- {({(u_{k})}_{\overline{x} x}, u_{l})}_{h} - λ_{k} {(u_{k}, u_{l})}_{h} = 0

Použijeme bu bun seki-bun s vědomím, že

{(u_{i})}_{0} = {(u_{i})}_{m} = 0

({(u_{k})}_{\overline{x}}, {(u_{l})}_{\overline{x}}] - λ_{k} {(u_{k}, u_{l})}_{h} = 0

a to samé samozřejmě platí i s prokozením

k \leftrightarrow l

({(u_{l})}_{\overline{x}}, {(u_{k})}_{\overline{x}}] - λ_{k} {(u_{l}, u_{k})}_{h} = 0

Odečtením dostaneme

(λ_{l} - λ_{k}) {(u_{k} - u_{l})}_{h} = 0

Takže pro

k \neq l

{(u_{k}, u_{l})}_{h} = 0

. Jelikož máme bázi, každou síťovou funkci

v \in ℋ_{h}

si můžeme vyjádřit jako lineární kombinaci:

v = \sum_{k = 1}^{m - 1} α_{k} u_{k}

Díky ortogonalitě máme jednoduché vyjádření pro normu:

{‖ v ‖}_{h}^{2} = {(\sum_{k = 1}^{m - 1} α_{k} u_{k}, \sum_{k = 1}^{m - 1} α_{k} u_{k})}_{h} = \sum_{k = 1}^{m - 1} \sum_{l = 1}^{m - 1} α_{k} α_{l} {(u_{k}, u_{l})}_{h} = \sum_{k = 1}^{m - 1} α_{k}^{2} {‖ u_{k} ‖}_{h}^{2}

To by ostatně plynulo přímo z Pythagorovy věty, tady jsme si vlastně rozepsali její důkaz. Analogicky si vyjádříme

{‖ v_{\overline{x}} ⟧}^{2} = \sum_{k = 1}^{m - 1} α_{k}^{2} ‖ u_{\overline{x}} ⟧ = \sum_{k = 1}^{m - 1} α_{k}^{2} λ_{k} {‖ u_{k} ‖}_{h}^{2}

Ze vzorce pro

λ_{k}

snadno odvodíme, že

λ_{1} < \dots < λ_{m - 1}

. Z toho plyne

λ_{1} \sum_{k = 1}^{m - 1} α_{k}^{2} \leq \sum_{k = 1}^{m - 1} α_{k}^{2} λ_{k} {‖ u_{k} ‖}_{h}^{2} \leq λ_{m - 1} \sum_{k = 1}^{m - 1} α_{k} {‖ u_{k} ‖}_{h}^{2}

neboli

λ_{1} {‖ v ‖}_{h}^{2} \leq ‖ v_{\overline{x}} ⟧ \leq λ_{m - 1} {‖ v ‖}_{h}^{2}

Snadno odhadneme

λ_{m - 1} = \frac{4}{h^{2}} {(\sin \frac{(m - 1) π}{2 m})}^{2} \leq \frac{4}{h^{2}}

, z čehož plyne první nerovnost. Druhá bude o něco složitější. Použijeme „metodu luku“. Nakreslíme si graf funkce

sin

na intervalu

(0, \frac{π}{4})

. Na tomto intervalu je funkce konkávní, z čehož plyne

\sin \frac{π}{2 m} \geq \frac{\frac{π}{2 m} \cdot \sin \frac{π}{4}}{\frac{π}{4}} = \frac{\sqrt{2}}{m}

λ_{1} = \frac{4}{h^{2}} {(\sin \frac{π}{2 m})}^{2} \geq \frac{8}{h^{2} m^{2}} = \frac{8}{{(b - a)}^{2}}

z čehož plyne druhá rovnost.

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_{\overline{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_{\overline{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)}_{\overline{x}})}_{x} + 𝒫_{h} ({(p y^{'})}^{'}), ψ_{0} = ψ_{m} = 0

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

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

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

- {(p z_{\overline{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_{\overline{x}})}_{x}, z)}_{h} + {(q z, z)}_{h} = - {(ψ, z)}_{h}

Použijeme bu bun seki-bun:

(p z_{\overline{x}}, z_{\overline{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_{\overline{x}}, z_{\overline{x}}] = \sum_{j = 1}^{m} h p_{j} {(z_{\overline{x}})}_{j}^{2} \geq c_{0} {‖ z_{\overline{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_{\overline{x}} ⟧}^{2} \leq {‖ ψ ‖}_{h} {‖ z ‖}_{h} \leq \frac{b - a}{\sqrt{8}} {‖ ψ ‖}_{h} ‖ z_{\overline{x}} ⟧

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

A z první Sobolevovy nerovnosti (společně s tou třetí) je

{‖ z ‖}_{0, h} \leq c_{L 1} ‖ z_{\overline{x}} ⟧ \leq c_{L 1} \frac{b - a}{\sqrt{8}} {‖ ψ ‖}_{h}

Z toho všeho a z Laxovy věty nějak plyne, ž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_{\overline{x} x} = f, x \in ω_{h}, - {(u_{x})}_{0} = β_{1} u_{0} + γ_{1}, {(u_{\overline{x}})}_{m} = β_{2} u_{m} + γ_{2}

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

\begin{aligned} - \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{aligned}

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{array}{cccc} - \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{array}), φ ≔ (\begin{array}{c} \frac{2 γ_{1}}{h} \\ f_{1} \\ ⋮ \\ f_{m - 1} \\ \frac{2 γ_{2}}{h} \end{array})

Uvažujme množinu síťových funkcí $ℋ_{h} ≔ {v : {\overline{ω}}_{h} \to ℝ}$ 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.

Věta.

𝐀

je samosdružený lineární operátor na

ℋ_{h}

, to jest

\forall u, v \in ℋ : [𝐀 u, v] = [u, 𝐀 v]

Věta. Existuje-li konstanta

C_{1} \in ℝ^{+}

taková, že

- β_{1} \geq C_{1} \land - β_{2} \geq C_{1}

, potom

𝐀

je pozitivně definitní, to jest

\forall u \in ℋ_{H} : [𝐀 u, u] \geq {(b - a)}^{2} \frac{1 + \frac{1}{C_{1} (b - a)}}{2} {⟦ u ⟧}^{2}

(ano, tohle vypadá naprosto přesně jako normální definice pozitivní definitnosti).

Věta. Schéma je stabilní, řešení diferenční rovnice konverguje k řešení diferenciální rovnice jako

𝒪 (\sqrt{h})

v normě

{⟦ \cdot ⟧}^{2}

a jako

𝒪 (h)

v normě

{‖ \cdot ‖}_{0, 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 : Ω \to ℝ, y ≔ y (x)$ .

Budeme předpokládat $f \in 𝒞 (Ω), g \in 𝒞 (\partial Ω), λ^{'}, q \in 𝒞 (Ω), q \geq 0, λ \geq C_{0} \in ℝ^{+}$ .

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

{\overline{ω}}_{h} ≔ {[i h_{1}, j h_{2}] \in Ω | 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_{{\overline{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_{{\overline{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_{{\overline{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_{{\overline{x}}_{1}})}_{x_{1}} - {(λ y_{{\overline{x}}_{2}})}_{x_{2}} + q u = f, x \in ω_{h}, u |_{γ_{h}} = g |_{γ_{h}}

Označíme-li ${\overline{\nabla}}_{h} u ≔ (u_{{\overline{x}}_{1}}, u_{{\overline{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 (λ {\overline{\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.

Věta (Greenova formule). Nechť

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

. Potom pro

λ : {\overline{ω}}_{h} \to ℝ, u, v \in ℋ_{h}

platí

{(\nabla_{h} \cdot (λ {\overline{\nabla}}_{h} u), v)}_{h} = - (λ {\overline{\nabla}}_{h} u, {\overline{\nabla}}_{h} v]

kde skalární součin definujeme jako

{(v, w)}_{h} ≔ \sum_{i = 1}^{m_{1} - 1} \sum_{j = 1}^{m_{2} - 1} h_{1} h_{2} u_{i, j} v_{i, j}

a ten druhý definujeme na konci důkazu.

Důkaz.

{(\nabla_{h} \cdot (λ {\overline{\nabla}}_{h} u), v)}_{h} = {({(λ u_{{\overline{x}}_{1}})}_{x_{1}} + {(λ u_{{\overline{x}}_{2}})}_{x_{2}}, v)}_{h} = {({(λ u_{{\overline{x}}_{1}})}_{x_{1}}, v)}_{h} + {({(λ u_{{\overline{x}}_{2}})}_{x_{2}}, v)}_{h}

Na oba výrazy použijeme jednorozměrnou Greenovu formuli:

\begin{aligned} {({(λ u_{{\overline{x}}_{1}})}_{x_{1}}, v)}_{h} & = \sum_{i = 1}^{m_{1} - 1} \sum_{j = 1}^{m_{2} - 1} h_{1} h_{2} {({(λ u_{{\overline{x}}_{1}})}_{x_{1}})}_{i, j} v_{i, j} \\ = \sum_{j = 1}^{m_{2} - 1} h_{2} \sum_{i = 1}^{m_{1} - 1} h_{1} {({(λ u_{{\overline{x}}_{1}})}_{x_{1}})}_{i, j} v_{i, j} \\ = \sum_{j = 1}^{m_{2} - 1} h_{2} \sum_{i = 1}^{m_{1}} h_{1} λ_{i, j} {(u_{{\overline{x}}_{1}})}_{i, j} {(v_{{\overline{x}}_{1}})}_{i, j} \\ = - \sum_{i = 1}^{m_{1}} \sum_{j = 1}^{m_{2} - 1} h_{1} h_{2} λ_{i, j} {(u_{{\overline{x}}_{1}})}_{i, j} {(v_{{\overline{x}}_{1}})}_{i, j} ≕ - (λ u_{{\overline{x}}_{1}}, v_{{\overline{x}}_{1}} ⌋ \end{aligned}

{({(λ u_{{\overline{x}}_{2}})}_{x_{2}}, v)}_{h} = \dots = - \sum_{i = 1}^{m_{1} - 1} \sum_{j = 1}^{m_{2}} h_{1} h_{2} λ_{i, j} {(u_{{\overline{x}}_{2}})}_{i, j} {(v_{{\overline{x}}_{2}})}_{i, j} ≕ - (λ u_{{\overline{x}}_{2}}, v_{{\overline{x}}_{2}} ⌉

Vždy, když mám pocit, že už Beneš nedokáže vymyslet divnější značení, tak mě přesvědčí o opaku. Teď to dáme dohromady:

{(\nabla_{h} \cdot (λ {\overline{\nabla}}_{h} u), v)}_{h} = - (λ u_{{\overline{x}}_{1}}, v_{{\overline{x}}_{1}} ⌋ - (λ u_{{\overline{x}}_{2}}, v_{{\overline{x}}_{2}} ⌉ ≕ - (λ {\overline{\nabla}}_{h} u, {\overline{\nabla}}_{h} v]

Lemma (druhá Sobolevova nerovnost). Nechť

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

. Potom existuje

C_{2,2} \in ℝ^{+}

takové, že pro všechna

u \in ℋ_{h}

{‖ u ‖}_{h} < C_{2,2} ‖ {\overline{\nabla}}_{h} u ⟧

, kde

‖ w ⟧ ≔ \sqrt{(w, w]}

Důkaz. Použijeme jednorozměrnou druhou Sobolevovu nerovnost v obou směrech.

{‖ u_{∙, j} ‖}_{h_{1}}^{2} = \sum_{i = 1}^{m_{1} - 1} h_{1} u_{i, j}^{2} \leq \frac{L_{1}^{2}}{4} \sum_{i = 1}^{m_{1}} h_{1} {(u_{{\overline{x}}_{1}})}_{i, j}^{2}

{‖ u_{i, ∙} ‖}_{h_{2}}^{2} = \sum_{j = 1}^{m_{2} - 1} h_{2} u_{i, j}^{2} \leq \frac{L_{2}^{2}}{4} \sum_{j = 1}^{m_{2}} h_{2} {(u_{{\overline{x}}_{2}})}_{i, j}^{2}

Vysčítáme první nerovnost přes všechny sloupce a druhou nerovnost přes všechny řádky. Obě vzniklé nerovnosti zprůměrujeme, čímž dostaneme

{‖ u ‖}_{h}^{2} \leq \frac{L_{1}^{2}}{8} (u_{{\overline{x}}_{1}}, u_{{\overline{x}}_{1}} ⌋ + \frac{L_{2}^{2}}{8} (u_{{\overline{x}}_{2}}, u_{{\overline{x}}_{2}} ⌉ \leq \frac{\max {L_{1}^{2}, L_{2}^{2}}}{8} ‖ {\overline{\nabla}}_{h} u ⟧

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 (λ {\overline{\nabla}}_{h} u) + q (𝒫_{h} y - u) = 0 = 0, 𝒫_{h} y |_{γ_{h}} = u |_{γ_{h}}

Použijeme chybu aproximace:

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

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

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

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

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

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

Aplikujeme Greenovu formuli:

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

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

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

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

Z energetické rovnosti a druhé Sobolevovy nerovnosti plyne

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

‖ {\overline{\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} ‖ {\overline{\nabla}}_{h} z ⟧ \leq \frac{c_{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:

{\overline{ω}}_{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} ≔ {\overline{ω}}_{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_{\overline{x} x})}_{j}^{k} + 𝒪 (h^{2})

\partial_{t} y |_{j}^{k} = {\begin{matrix} {(y_{t})}_{j}^{k} + 𝒪 (τ) = \frac{y_{j}^{k + 1} - y_{j}^{k}}{τ} + 𝒪 (τ) \\ {(y_{\overline{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{matrix}

Explicitní schéma

Použijeme dopřednou diferenci.

\frac{\hat{u} - u}{τ} = D u_{\overline{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_{\overline{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{array}{cccc} 1 - \frac{2 τ D}{h^{2}} & \frac{τ d}{h^{2}} & 0 & \dots \\ \frac{τ d}{h^{2}} & ⋱ & ⋱ & ⋱ \\ 0 & ⋱ & ⋱ & ⋱ \end{array}]

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:

{\overline{λ}}_{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}}_{\overline{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}}_{\overline{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{array}{cccc} 1 + \frac{2 τ D}{h^{2}} & - \frac{τ d}{h^{2}} & 0 & \dots \\ - \frac{τ d}{h^{2}} & ⋱ & ⋱ & ⋱ \\ 0 & ⋱ & ⋱ & ⋱ \end{array}]

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

U^{k + 1} = 𝐁^{- 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}}_{\overline{x} x} + u_{\overline{x} x}}{2} + \frac{\hat{f} + f}{2}

Chyba aproximace je opět $𝒪 (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_{\overline{x} x})}_{j}^{k + 1} + {(u_{\overline{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_{\overline{x} x})}^{k + 1} + {(y_{\overline{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_{\overline{x} x})}^{k + 1} + {(y_{\overline{x} x})}^{k}) + \frac{D}{2} ({(u_{\overline{x} x})}_{j}^{k + 1} + {(u_{\overline{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}

z ≔ 𝒫_{h} y - u

???????????????????????????????????????????

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

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

Bu bun seki-bun asi?

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

TBD

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_{\overline{x} x} + 𝒪 (h^{2})

\dot{u} = D u_{\overline{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:

{\dot{u}}_{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 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í

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

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

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