Výpisky z Numerické matematiky 1

Stacionární iterativní metody

Obecný předpis: ${\vec{x}}^{(k + 1)} = 𝐁^{(k)} {\vec{x}}^{(k)} + {\vec{c}}^{(k)}$ , přičemž přesné řešení ${\vec{x}}^{*}$ by mělo být pevný bod pro libovolné $k$ .

Věta. Iterativní metoda konverguje z libovolného počátečního bodu ke svému pevnému bodu právě tehdy, pokud

\lim_{k \to \infty} \prod_{0}^{i = k} 𝐁^{(i)} = 0

Důkaz.

{\vec{x}}^{(k + 1)} - {\vec{x}}^{*} = 𝐁^{(k)} {\vec{x}}^{k} + {\vec{c}}^{(k)} - 𝐁^{*} {\vec{x}}^{k} - {\vec{c}}^{(k)} = 𝐁^{(k)} ({\vec{x}}^{(k)} - {\vec{x}}^{*})

Tím dostáváme jakýsi rekurentní vztah. Rozvinutím máme

{\vec{x}}^{(k + 1)} - {\vec{x}}^{*} = \prod_{0}^{i = k} 𝐁^{(i)} ({\vec{x}}^{(0)} - {\vec{x}}^{*})

což jsme chtěli dokázat.

Díky tomuto jsou iterativní metody samoopravující — i když ve výpočtu uděláme chybu, v dalších krocích se budeme opět přibližovat

Definice. Iterativní metoda je stacionární, pokud koeficienty

𝐁, \vec{c}

nezávisí na čísle kroku.

Důsledek. Stacionární iterativní metoda konverguje z libovolného počátečního bodu ke svému pevnému bodu právě tehdy, pokud

\lim_{k \to \infty} 𝐁^{k} = 0

Důsledek. Stacionární iterativní metoda konverguje z libovolného počátečního bodu ke svému pevnému bodu právě tehdy, pokud

ρ (𝐁) < 1

Důsledek. Stacionární iterativní metoda konverguje z libovolného počátečního bodu ke svému pevnému bodu, pokud existuje taková norma, že

‖ 𝐁 ‖ < 1

Věta (aposteriorní odhad pomocí residua). Pro stacionární iterativní metodu k řešení soustavy lineárních rovnic platí následující odhad chyby:

‖ {\vec{x}}^{(k)} - {\vec{x}}^{*} ‖ \leq ‖ 𝐀^{- 1} ‖ ‖ 𝐀 {\vec{x}}^{(k)} - \vec{b} ‖ = ‖ 𝐀^{- 1} ‖ ‖ {\vec{r}}^{(k)} ‖

kde maticová norma je souhlasná s vektorovou normou.

Důkaz.

{\vec{x}}^{(k)} - {\vec{x}}^{*} = 𝐀^{- 1} (𝐀 {\vec{x}}^{(k)} - \vec{b}) = 𝐀^{- 1} {\vec{r}}^{(k)}

Věta (aposteriorní odhad pomocí změny). Pro jakoukoli stacionární iterativní metodu platí následující odhad chyby:

‖ {\vec{x}}^{(k)} - {\vec{x}}^{*} ‖ \leq ‖ {(𝐈 - 𝐁)}^{- 1} ‖ ‖ {\vec{x}}^{(k - 1)} - {\vec{x}}^{(k)} ‖

kde maticová norma je souhlasná s vektorovou normou.

Důkaz. Nejprve dokážeme, že

𝐈 - 𝐁

je regulární. Rozepíšeme si matici

𝐁

podle Schurovy věty jako

𝐁 ≕ 𝐔^{*} 𝐑 𝐔

. Víme, že pokud

{\vec{x}}^{(k)}

konverguje, potom

ρ (𝐁) < 1

, tedy

𝐑

má na diagonále čísla absolutně menší než

1

. V takovém případě matice

𝐈 - 𝐑

je trojúhelníková s nenulovými prvky na diagonále, a tedy regulární, tudíž i

𝐈 - 𝐁 = 𝐔^{*} (𝐈 - 𝐑) 𝐔

je regulární. Nyní můžeme psát:

{\vec{x}}^{*} = 𝐁 {\vec{x}}^{*} + \vec{c} ∴ (𝐈 - 𝐁) {\vec{x}}^{*} = \vec{c} ∴ x^{*} = {(𝐈 - 𝐁)}^{- 1} \vec{c}

\begin{aligned} {\vec{x}}^{(k)} - {\vec{x}}^{*} & = {\vec{x}}^{(k)} - {(𝐈 - 𝐁)}^{- 1} \vec{c} \\ = {(𝐈 - 𝐁)}^{- 1} ((𝐈 - 𝐁) {\vec{x}}^{(k)} - \vec{c}) \\ = {(𝐈 - 𝐁)}^{- 1} ({\vec{x}}^{(k)} - 𝐁 {\vec{x}}^{(k)} - \vec{c}) \\ = {(𝐈 - 𝐁)}^{- 1} ({\vec{x}}^{(k)} - {\vec{x}}^{(k + 1)}) \end{aligned}

z čehož plyne tvrzení věty.

Věta (apriorní odhad). Pro stacionární iterativní metodu platí následující odhad chyby:

‖ {\vec{x}}^{(k)} - {\vec{x}}^{*} ‖ \leq {‖ 𝐁 ‖}^{k} (‖ {\vec{x}}^{(0)} ‖ + \frac{‖ \vec{c} ‖}{1 - ‖ 𝐁 ‖})

kde maticová norma je souhlasná s vektorovou normou.

Důkaz. Z důkazu předchozí věty máme:

{\vec{x}}^{*} = {(𝐈 - 𝐁)}^{- 1} \vec{c} = \sum_{i = 0}^{\infty} 𝐁^{i} \vec{c}

Nyní si vyjádříme:

{\vec{x}}^{(k)} = 𝐁 {\vec{x}}^{(k - 1)} + \vec{c} = 𝐁 (𝐁 {\vec{x}}^{(k - 2)} + \vec{c}) + \vec{c} = \dots = 𝐁^{k} {\vec{x}}^{(0)} + \sum_{i = 0}^{k - 1} 𝐁^{i} \vec{c}

{\vec{x}}^{(k)} - {\vec{x}}^{*} = 𝐁^{k} {\vec{x}}^{(0)} + \sum_{i = 0}^{k - 1} 𝐁^{i} \vec{c} - \sum_{i = 0}^{\infty} 𝐁^{i} \vec{c} = 𝐁^{k} {\vec{x}}^{(0)} - \sum_{i = k}^{\infty} 𝐁^{i} \vec{c} = 𝐁^{k} ({\vec{x}}^{(0)} - \sum_{i = 0}^{\infty} 𝐁^{i} \vec{c})

Použitím několika trojúhelníkových nerovností a sečtením řady dostaneme tvrzení věty.

Definice. Metoda postupných aproximací je metoda pro řešení soustavy lineárních rovnic založená na odečítání rezidua:

{\vec{x}}^{(k + 1)} ≔ {\vec{x}}^{(k)} - {\vec{r}}^{(k)} = (𝐈 - 𝐀) {\vec{x}}^{(k)} + \vec{b}

Důsledek. Metoda postupných aproximací konverguje z libovolného počátečního bodu ke svému pevnému bodu právě tehdy, pokud

ρ (𝐈 - 𝐀) < 1

Důsledek. Metoda postupných aproximací konverguje z libovolného počátečního bodu ke svému pevnému bodu, pokud existuje taková norma, že

‖ 𝐈 - 𝐀 ‖ < 1

Věta. Nechť

p

je polynom a

𝐀

čtvercová matice. Potom

(\forall λ \in σ (𝐀)) (p (λ) \in σ (p (𝐀)))

Důkaz.

𝐀 \vec{x} = λ \vec{x} ⟹ p (𝐀) \vec{x} = \sum_{i = 0}^{s} p_{i} 𝐀^{i} \vec{x} = \sum_{i = 0}^{s} p_{i} λ^{i} \vec{x} = p (λ) \vec{x}

Věta. Je-li

𝐀

hermitovská, potom metoda postupných aproximací konverguje právě tehdy, pokud

𝟎 < 𝐀 < 2 𝐈

Důkaz.

\begin{aligned} metoda konverguje & ⟺ ρ (𝐁) < 1 \\ ⟺ ρ (𝐈 - 𝐀) < 1 \\ ⟺ σ (𝐈 - 𝐀) \subseteq (- 1,1) \\ ⟺ (\forall λ \in σ (𝐀)) (1 - λ \in (- 1,1)) \\ ⟺ (\forall λ \in σ (𝐀)) (λ \in (0,2)) \\ ⟺ 𝐀 > 𝟎 \land (\forall λ \in σ (𝐀)) (λ < 2) \\ ⟺ 𝐀 > 𝟎 \land (\forall λ \in σ (𝐀)) (2 - λ > 0) \\ ⟺ 𝐀 > 𝟎 \land 2 𝐈 - 𝐀 > 0 \end{aligned}

To znamená, že matic, pro které konverguje, není moc. S tím nám pomůže předpodmíněni — vynásobíme obě strany rovnice vhodnou regulární maticí $𝐇$ .

Věta. Je-li

𝐀

hermitovská a pozitivně definitní, potom předpodmíněná metoda postupných aproximací konverguje právě tehdy, pokud

𝟎 < 𝐀 < 𝐖 + 𝐖^{*}

, kde

𝐖 ≔ 𝐇^{- 1}

. V takovém případě je konvergence monotónní vzhledem k energetické normě

{‖ \cdot ‖}_{𝐀}

Důkaz.

\begin{aligned} {‖ 𝐁 ‖}_{𝐀} & = {‖ 𝐀^{\frac{1}{2}} {𝐁 𝐀}^{- \frac{1}{2}} ‖}_{2} \\ = {‖ 𝐀^{\frac{1}{2}} (𝐈 - 𝐇 𝐀) 𝐀^{- \frac{1}{2}} ‖}_{2} \\ = {‖ 𝐈 - 𝐀^{\frac{1}{2}} 𝐇 𝐀^{\frac{1}{2}} ‖}_{2} \\ ≕ {‖ \hat{𝐁} ‖}_{2} \\ = \sqrt{ρ ({\hat{𝐁}}^{*} \hat{𝐁})} \end{aligned}

Obecně pro libovolnou matici

𝐗

platí, že

𝐗^{*} 𝐗

je hermitovská a pozitivně definitní (důkaz triviální). Takže stačí dokázat, že všechna vlastní čísla

{\hat{𝐁}}^{*} \hat{𝐁}

jsou menší než

1

\begin{aligned} {\hat{𝐁}}^{*} \hat{𝐁} & = {(𝐈 - 𝐀^{\frac{1}{2}} 𝐇 𝐀^{\frac{1}{2}})}^{*} (𝐈 - 𝐀^{\frac{1}{2}} 𝐇 𝐀^{\frac{1}{2}}) \\ = (𝐈 - 𝐀^{\frac{1}{2}} 𝐇^{*} 𝐀^{\frac{1}{2}}) (𝐈 - 𝐀^{\frac{1}{2}} 𝐇 𝐀^{\frac{1}{2}}) \\ = 𝐈 - 𝐀^{\frac{1}{2}} (𝐇^{*} + 𝐇) 𝐀^{\frac{1}{2}} + 𝐀^{\frac{1}{2}} 𝐇^{*} 𝐀 𝐇 𝐀^{\frac{1}{2}} \\ = 𝐈 - 𝐀^{\frac{1}{2}} 𝐇^{*} (𝐇^{- 1} + {𝐇^{*}}^{- 1} - 𝐀) 𝐇 𝐀^{\frac{1}{2}} \\ = 𝐈 - 𝐀^{\frac{1}{2}} 𝐇^{*} (𝐖 + 𝐖^{*} - 𝐀) 𝐇 𝐀^{\frac{1}{2}} \end{aligned}

Pokud platí předpoklad, matice v závorce je pozitivně definitní. Dá se dokázat, že pro libovolné matice

𝐗, 𝐘

platí

𝐗 > 𝟎 ⟹ 𝐘^{*} 𝐗 𝐘 > 𝟎

, z čehož plyne, že vlastní čísla původního výrazu jsou menší než

1

. Tudíž energetická norma

𝐁

je menší než

1

a metoda tedy konverguje. Dokazovat zbytek věty se Oberhuberovi nechtělo.

Dá se snadno ukázat, že by bylo hezké, kdyby $𝐇$ byla co nejblíž $𝐀^{- 1}$ . Ovšem počítat inverzi samozřejmě nepřipadá v úvahu, takže se ji budeme snažit aproximovat. Dá se použít například takzvaný neúplný LU rozklad (ale v praxi se zrovna pro tohle tolik nepoužívá).

Definice. Nechť

θ \in ℝ

je relaxační parametr. Richardsonovy iterace jsou iterační metoda pro řešení soustavy lineárnich rovnic daná předpisem

{\vec{x}}^{(k + 1)} ≔ {\vec{x}}^{(k)} - θ {\vec{r}}^{(k)} = {\vec{x}}^{(k)} - θ (𝐀 {\vec{x}}^{(k)} - \vec{b})

Poznámka. Jde o metodu postupných aproximací s předpodmíněním

θ 𝐈

Důsledek. Je-li

𝐀

hermitovská a pozitivně definitní, potom Richardsonova metoda konverguje právě tehdy, pokud

𝐀 < \frac{2}{θ} 𝐈

Poznámka. Chceme tedy zvolit

θ < \frac{2}{ρ (𝐀)}

Jacobiho metoda: v každém kroku vezmeme $i$ -tou rovnici, vyjádříme z ní $x_{i}$ a za ostatní dosadíme z předchozího kroku. Jak to vyjádřit formálně? Zapíšeme si $𝐀 ≕ 𝐃 - 𝐋 - 𝐑$ , kde trojúhelníkové matice mají nuly na diagonále a provedeme jednoduché úpravy, z čehož vznikne vzorec

{\vec{x}}^{(k + 1)} ≔ 𝐃^{- 1} (𝐋 + 𝐑) {\vec{x}}^{(k)} + 𝐃^{- 1} \vec{b} = (𝐈 - 𝐃^{- 1} 𝐀) {\vec{x}}^{(k)} + 𝐃^{- 1} \vec{b}

To je vlastně metoda postupných aproximací s podmíněním $𝐃^{- 1}$ , kterému se zove Jacobiho podmínění.

Věta. Jacobiho metoda konverguje právě tehdy, pokud

ρ (𝐃^{- 1} (𝐋 + 𝐑)) < 1

Definice. Matice s převládající diagonálou je matice, kde v každém řádku je absolutní hodnota diagonálního prvku ostře větší než součet absolutních hodnot ostatních prvků.

Věta. Je-li

𝐀

matice s převládající diagonálou, potom Jacobiho metoda konverguje.

Důkaz. Dokážeme, že

{‖ 𝐁 ‖}_{\infty} < 1

. Vyjádříme si, jak vypadá

𝐁

𝐁 = 𝐃^{- 1} (𝐋 + 𝐑)

B_{i, j} = {\begin{matrix} 0, & i = j \\ - \frac{A_{i, j}}{A_{i, i}}, & i \neq j \end{matrix}

{‖ 𝐁 ‖}_{\infty} = \max_{i \in \hat{n}} \sum_{j = 1}^{n} | B_{i, j} | = \max_{i \in \hat{n}} \sum_{\begin{array}{c} j = 1 \\ j \neq i \end{array}}^{n} \frac{| A_{i, j} |}{| A_{i, i} |} < 1

Věta. Je-li

𝐀

hermitovská a pozitivně definitní, potom Jacobiho metoda konverguje právě tehdy, pokud

𝐀 < 2 𝐃

. V takovém případě je konvergence monotónní vzhledem k energetické normě

{‖ \cdot ‖}_{𝐀}

Důkaz. Použijeme obecnou větu o konvergenci předpodmíněných metod, kde

𝐖 = 𝐃

. Jelikož

𝐀

je hermitovská, máme

𝐖 + 𝐖^{*} = 𝐃 + 𝐃^{*} = 2 𝐃

Jacobiho metodu můžeme vylepšit tím, že při výpočtu budeme rovnou používat již napočítané složky, čímž získáme Gaussovu-Seidelovu metodu. Jde o metodu postupných aproximací s předpodmíněním ${(𝐃 - 𝐋)}^{- 1}$ .

Věta. Je-li

𝐀

matice s převládající diagonálou, potom Gaussova-Seidelova metoda konverguje.

Důkaz.

𝐁 = (𝐃 - 𝐋^{- 1}) 𝐑

{‖ 𝐁 ‖}_{\infty} = \max {‖ 𝐁 \vec{x} ‖ | ‖ \vec{x} ‖ = 1}

Nechť

\vec{u}

je vektor, pro který je dosaženo maxima, a

\vec{v} ≔ 𝐁 \vec{u}

, tedy

{‖ 𝐁 ‖}_{\infty} = {‖ \vec{v} ‖}_{\infty}

. Nechť

s

je index absolutně nejvyšší složky

\vec{v}

, tedy

{‖ 𝐁 ‖}_{\infty} = | v_{s} |

. Podle definice

𝐁

máme

(𝐃 - 𝐋) \vec{v} = 𝐑 \vec{u}

, neboli

(𝐃 - 𝐋) \vec{v} - 𝐑 \vec{u} = \vec{0}

. To je nějaká soustava rovnic, podíváme se na její

s

-tou složku. Všimněme si, že ji umíme vyjádřit pomocí matice

𝐀

\sum_{i = 1}^{s} A_{s, i} v_{i} + \sum_{i = s + 1}^{n} A_{s, i} u_{i} = 0

Z toho vyjárďíme

v_{i}

a dáme ho do absolutní hodnoty:

\begin{aligned} ‖ 𝐁 ‖ & = | v_{i} | \\ = | - \frac{1}{A_{s, s}} (\sum_{i = 1}^{s - 1} A_{s, i} v_{i} + \sum_{i = s + 1}^{n} A_{s, i} u_{i} = 0) | \\ \leq \sum_{i = 1}^{s - 1} \frac{| A_{s, i} |}{| A_{s, s} |} | v_{i} | + \sum_{i = s + 1}^{n} \frac{| A_{s, i} |}{| A_{s, s} |} | u_{i} | \\ \leq | v_{s} | \sum_{i = 1}^{s - 1} \frac{| A_{s, i} |}{| A_{s, s} |} + \sum_{i = s + 1}^{n} \frac{| A_{s, i} |}{| A_{s, s} |} & j e l i k o \overset{ˇ}{z} {‖ u ‖}_{\infty} = 1, {‖ v ‖}_{\infty} = v_{s} \\ = a | v_{s} | + b, a + b = 1 & j e l i k o \overset{ˇ}{z} m a t i c e m \overset{ˊ}{a} p \overset{ˇ}{r} e v l \overset{ˊ}{a} d a j \overset{ˊ}{ı} c \overset{ˊ}{ı} d i a g o n \overset{ˊ}{a} l u \\ = \frac{b}{1 - a} < \frac{b}{b} = 1 \end{aligned}

Věta. Je-li

𝐀

hermitovská a pozitivně definitní, potom Gaussova-Seidelova metoda konverguje vždy. V takovém případě je konvergence monotónní vzhledem k energetické normě

{‖ \cdot ‖}_{𝐀}

Důkaz. Použijeme obecnou větu o konvergenci předpodmíněných metod, kde

𝐖 = 𝐃 - 𝐋

. Z hermitovskosti

𝐀

máme

𝐃^{*} = 𝐃, 𝐋^{*} = 𝐑, 𝐑^{*} = 𝐑

. Tudíž máme:

𝐖 + 𝐖^{*} = 𝐃 + 𝐃 + 𝐋 - 𝐑 = 𝐃 + 𝐀

Jelikož

𝐃

je pozitivně definitní, máme

𝟎 < 𝐀 < 𝐖 + 𝐖^{*} .

Poznámka. Gaussova-Seidelova metoda je obecně lepší než Jacobiho metoda, ale není jednoznačně lepší.

Pokud ke Gaussově-Seidelově metodě přidáme relaxační parametr, dostaneme super-relaxační metodu. Ta je vhodná pro lineární soustavy pocházející z metody konečných diferencí pro řešení parabolických nebo eliptických parciálních diferenciálních rovnic. Abychom pro ni dostali předpis, musíme vztah pro Gaussovu-Seidelovu metodu vyjádřit jako $x_{i}^{(k + 1)} = x_{i}^{(k)} + Δ x_{i}^{(k)}$ . Poté zavedeme relaxační parametr $ω \in ℝ$ a položíme $x_{i}^{(k + 1)} ≔ x_{i}^{(k)} + ω Δ x_{i}^{(k)}$ . Je to metoda postupných aproximací s podmíněním $ω {(𝐃 - ω 𝐋)}^{- 1}$ .

Věta. Pro super-relaxační metodu platí

ρ (𝐁) \geq | ω - 1 |

, tedy metoda nikdy nekonverguje pro

ω \notin (0,2)

Důkaz. Využijeme toho, že determinant matice se rovná součinu vlastních čísel (dá se dokázat pomocí Jordanova tvaru). Matice

𝐁

se dá vyjádřít ve tvaru:

𝐁 = {(𝐃 - ω 𝐋)}^{- 1} ((1 - ω) 𝐃 + ω 𝐑)

\det (𝐁) = \frac{\det ((1 - ω) 𝐃 + ω 𝐑)}{\det (𝐃 - ω 𝐋)}

Obě matice, jejichž determinant počítáme, jsou trojúhelníkové, takže stačí vynásobit prvky na diagonále. Máme:

\det (𝐁) = \prod_{i = 1}^{n} \frac{(1 - ω) D_{i, i}}{D_{i, i}} = {(1 - ω)}^{n}

Aby se součin nějakých čísel mohl rovnat něčemu, tak alespoň jedno musí být absolutně větší nebo rovné.

Věta. Je-li

𝐀

matice s převládající diagonálou, potom super-relaxační metoda konverguje pro libovolné

ω \in (0,1]

Důkaz. Neznámý.

Věta. Je-li

𝐀

hermitovská a pozitivně definitní, potom pro

ω \in (0,2)

super-relaxační metoda konverguje. V takovém případě je konvergence monotónní vzhledem k energetické normě

{‖ \cdot ‖}_{𝐀}

Důkaz. Opět použijeme svou oblíbenou větu, kde máme

𝐖 = ω^{- 1} (𝐃 - ω 𝐋)

𝐖 + 𝐖^{*} = ω^{- 1} (𝐃 - ω 𝐋) + ω^{- 1} (𝐃 - ω 𝐑) = \frac{2}{ω} 𝐃 - 𝐋 - 𝐑 = (\frac{2}{ω} - 1) 𝐃 + 𝐀

Chceme tedy, aby bylo

(\frac{2}{ω} - 1) 𝐃 > 𝟎

, což je splněno pro

ω \in (0,2)

Budeme pro tzv. dvoucyklické shodně uspořádané matice (ty se často vyskytují při řešení diferenciálních rovnic) zkoumat optimální hodnotu $ω$ .

Definice. Čtvercová matice

𝐂

je slabě cyklická s indexem 2, pokud existuje taková permutace

𝐏

, že

{𝐏 𝐂 𝐏}^{T} = (\begin{array}{cc} 𝟎 & 𝐌_{1} \\ 𝐌_{2} & 𝟎 \end{array})

Definice. Matice

𝐀 = 𝐃 - 𝐋 - 𝐑

se nazývá dvoucyklická, pokud její Jacobiho matice

𝐃^{- 1} (𝐋 + 𝐑)

je slabě cyklická s indexem 2. Dvoucyklická matice je shodně uspořádaná, pokud vlastní číslo matice

α 𝐃^{- 1} 𝐋 + α^{- 1} 𝐃^{- 1} 𝐑

je stejné pro všechna

α \in ℝ_{\neq 0}

Věta. Nechť

𝐀

je dvoucyklická shodně uspořádaná matice,

ω \in ℝ_{\neq 0}

λ \neq 0

je vlastní číslo matice

𝐁_{ω}

ze super-relaxační metody. Nechť pro nějaké

μ

platí

{(λ + μ - 1)}^{2} = ω^{2} μ^{2} λ

. Potom

μ

je vlastní číslo matice

𝐁_{J}

z Jacobiho metody, a vopáčně. Navíc

ρ (𝐁_{ω})

je nejmenší, a super-relaxační metoda tedy konverguje nejrychleji, pro

ω_{opt} ≔ \frac{2}{1 + \sqrt{1 - ρ {(𝐁_{J})}^{2}}}

Tohle je ale spíš teoretický výsledek, protože je většinou rychlejší dát si pár iterací navíc než počítat nějaký spektrální poloměr. Často se dělá to, že se parametrem $ω$ pohybuje a hledá se, pro jakou hodnotu to konverguje nejrychleji. Ukazuje se, že se vyplatí ho hledat přesně (až na dvě desetinná místa).

Mocninná metoda

Mocninná metoda řeší částečný problém vlastních čísel — nenajde všechny, ale jen to absolutně největší, což je často jediné, co nás zajímá.

Základní myšlenka je taková, že vezmeme nějaký vektor $\vec{x}$ , budeme počítat posloupnost $𝐀^{k} \vec{x}$ (Krylovovu posloupnost) a sledovat, jak rychle roste. Je potřeba ho zvolit tak, aby měl nenulový průmět ve směru největšího vlastního vektoru. Pokud vektor v každém kroku posloupnosti znormalizuujeme (nejčastěji maximovou normou), aby posloupnost neutekla do nekonečna, dostáváme právě mocninnou metodu. Za vhodných podmínek potom bude konvergovat k vlastnímu vektoru a jeho největší složka $ρ_{k}$ k největšímu vlastnímu číslu. Pokud ovšem počáteční vektor zvolíme blbě, můžeme dostat jiné vlastní číslo.

Věta. Nechť matice

𝐀

má jedno absolutně největší číslo

λ

se stejnou algebraickou i geometrickou násobností

r

. Nechť její Jordanův tvar je

𝐗^{- 1} 𝐉 𝐗

, kde

𝐉

začíná

r

řádky s

λ

. Potom pro libovolný počáteční odhad

{\vec{x}}^{(0)}

takový, že

{(𝐗 {\vec{x}}^{(0)})}_{1, \dots, r} \neq \vec{0}

, mocninná metoda najde vlastní číslo

λ

jako limitu

ρ_{k}

a nějaký příslušný vlastní vektor jako limitu

{\vec{x}}^{(k)}

Poznámka. Pokud jsou absolutně největší vlastní čísla dvě (s opačným znaménkem)

\pm λ

, konverguje k

| λ |

posloupnost

\sqrt{ρ_{2 k} ρ_{2 k + 1}}

a k vlastnímu vektoru příslušejícímu k

\pm λ

posloupnost

𝐀 {\vec{x}}^{(2 k)} \pm λ {\vec{x}}^{(2 k)}

Důkaz. Pro zjednodušení nebudeme vektor dělit maximem, ale první složkou. Předpis algoritmu tedy bude:

{\vec{x}}^{(k + 1)} ≔ \frac{𝐀 {\vec{x}}^{(k)}}{ρ_{k} ≔ {\vec{e}}_{1}^{T} 𝐀 {\vec{x}}^{(k)}}

\begin{aligned} ρ_{k} & = {\vec{e}}_{1}^{T} 𝐀 {\vec{x}}^{(k)} \\ = \frac{{\vec{e}}_{1}^{T} 𝐀 {\vec{x}}^{(k)}}{{\vec{e}}_{1}^{T} {\vec{x}}^{(k)}} \\ = \frac{{\vec{e}}_{1}^{T} 𝐀 (\prod_{i = 1}^{k} ρ_{i}^{- 1}) 𝐀^{k} {\vec{x}}^{(0)}}{{\vec{e}}_{1}^{T} (\prod_{i = 1}^{k} ρ_{i}^{- 1}) 𝐀^{k} {\vec{x}}^{(0)}} \\ = \frac{{\vec{e}}_{1}^{T} 𝐀^{k + 1} {\vec{x}}^{(0)}}{{\vec{e}}_{1}^{T} 𝐀^{k} {\vec{x}}^{(0)}} \\ = \frac{{\vec{e}}_{1}^{T} 𝐗^{- 1} 𝐉^{k + 1} 𝐗 {\vec{x}}^{(0)}}{{\vec{e}}_{1}^{T} 𝐗^{- 1} 𝐉^{k} 𝐗 {\vec{x}}^{(0)}} \\ ≕ \frac{{\vec{e}}_{1}^{T} 𝐗^{- 1} λ^{k + 1} {\tilde{𝐉}}^{k + 1} 𝐗 {\vec{x}}^{(0)}}{{\vec{e}}_{1}^{T} 𝐗^{- 1} λ^{k} {\tilde{𝐉}}^{k} 𝐗 {\vec{x}}^{(0)}} \end{aligned}

Přeškálovaná Jordanova matice bude vypadat takto:

\tilde{𝐉} = (\begin{array}{c} 1 \\ ⋱ \\ 1 \\ \frac{𝐉_{2}}{λ} \\ ⋱ \\ \frac{𝐉_{s}}{λ} \end{array})

kde každý blok má zřejmě spektrální poloměr menší než

1

, tudiž jeho mocnina bude konvergovat k nulové matici. Ve zlomku se tak všechno pokrátí a máme

\lim_{k \to \infty} ρ_{k} = λ

. Když podobným způsobem upravíme výraz pro

{\vec{x}}^{(k)}

, zjistíme, že konverguje k nějakému

\vec{x}

. Zlimitíme-li nyní předpis, dostáváme

\vec{x} = \frac{𝐀 \vec{x}}{λ}

tedy jde o vlastní vektor.

Rychlost konvergence závisí na poměru hledaného vlastního čísla a ostatních vlastních čísel. Metodu můžeme urychlit tím, že spektrum vhodně posuneme (od $𝐀$ odečteme nějaké $λ^{*} 𝐈$ ).

Pokud chceme navopák najít nejmenší vlastní číslo, stačí použít metodu na matici $𝐀^{- 1}$ . Tu ale nechceme počítat, takže místo toho budeme vždy řešit soustavu $𝐀 {\vec{x}}^{(k + 1)} = {\vec{x}}^{(k)}$ nějakou iterativní metodou. To půjde rychle, protože pro malá $k$ nám stačí malá přesnost a pro velká $k$ se vektor mění málo, takže máme dobrý počáteční odhad.

Pokud chceme vlastních čísel najít víc, můžeme použít redukční metodu, která z matice „odstraní“ jedno vlastní číslo. K výpočtu kompletního spektra se ale nehodí, protože způsobuje ztrátu přesnosti. Nechť má matice vlastní vektor $\vec{x}$ . Pomocí matice přechodu $𝐏$ ji převedeme do báze $\vec{x}, {\vec{e}}_{2}, \dots, {\vec{e}}_{n}$ :

𝐏 = (\begin{array}{c} x_{1} \\ x_{2} & 1 \\ ⋮ & ⋱ \\ x_{n} & 1 \end{array})

𝐏^{- 1} = (\begin{array}{c} \frac{1}{x_{1}} \\ - \frac{x_{2}}{x_{1}} & 1 \\ ⋮ & ⋱ \\ - \frac{x_{n}}{x_{1}} & 1 \end{array})

𝐏^{- 1} 𝐀 𝐏 = (\begin{array}{cc} λ & {\vec{q}}^{T} \\ \vec{0} & 𝐁 \end{array})

Matice $𝐁$ bude mít stejná vlastní čísla jako $𝐀$ , ale bez $λ$ , takže můžeme mocninnou metodou najít nějaké další vlastní číslo $μ$ s vlastním vektorem $\vec{z}$ . Ale jak příslušný vlastní vektor převedeme na vlastní vektor $𝐀$ ? Nejprve mu dopočteme první složku

z_{1} ≔ \frac{{\vec{q}}^{T} \vec{z}}{μ - λ}

Pokud $λ = μ$ , vlastní číslo je vícenásobné, takže ${\vec{z}}_{1}$ můžu volit libovolně. Ze $\vec{z}$ dopočtu vlastní vektor $\vec{y}$ matice $𝐀$ jako $\vec{y} = 𝐏 \vec{z}$ .

Trojúhelníková metoda, LR algoritmus

Co když chceme najít všechna vlastní čísla? Myšlenka trojúhelníkové metody je taková, že použijeme mocninnou metodu na více vektorů, ale budeme zajišťovat, aby zůstaly lineárně nezávislé. Budeme konstruovat posloupnosti dolních a horních trojúhelníkových matic. $𝐋^{(0)}$ zvolíme libovolně a následně budeme brát $𝐋^{(k + 1)} 𝐑^{(k + 1)} ≔ {𝐀 𝐋}^{(k)}$ . Pokud tyto posloupnosti konvergují, pak platí $𝐀 𝐋 = 𝐋 𝐑$ , tedy matice $𝐑$ je podobná matici $𝐀$ a na její diagonále najdeme vlastní čísla. Vlastní vektory $𝐑$ potom najdeme řešením soustavy $(𝐑 - λ 𝐈) \vec{x} = \vec{0}$ , ty potom přenásobíme $𝐋$ a dostaneme vlastní vektory $𝐀$ .

$𝐋^{(0)}$ ani nemusí být dolní trojúhelníková, ale ${𝐀 𝐋}^{(0)}$ musí být silně regulární, abychom mohli provést LU rozklad. Jelikož můžeme začít z prakticky libovolné matice, metoda má samoopravující schopnosti. Zároveň potřebujeme, aby ${𝐀 𝐋}^{(k)}$ bylo silně regulární pro každé $k$ .

Věta. Je-li

𝐀

silně regulární, potom všechny matice na nějakém okolí jsou silně regulární.

Důkaz. Vzorec pro LU rozklad je spojitý, takže pro každou matici na nějakém okolí také existuje.

Věta. Nechť

𝐀 = 𝐈 + 𝐄 ≕ 𝐋 𝐑

. Potom

\lim_{‖ 𝐄 ‖ \to 0} 𝐋 = \lim_{‖ 𝐄 ‖ \to 0} 𝐑 = 𝐈

Důkaz. Opět plyne ze spojitosti.

Lemma. Mějme matice z trojúhelníkové metody. Existuje-li LU rozklad

𝐀^{k} 𝐋^{(0)} ≕ ℒ^{(k)} ℛ^{(k)}

, potom platí

ℒ^{(k)} = 𝐋^{(k)}

ℛ^{(k)} = \prod_{1}^{i = k} 𝐑^{(i)}

Důkaz.

𝐀^{k} 𝐋^{(0)} = 𝐀^{k - 1} 𝐋^{(1)} 𝐑^{(1)} = 𝐀^{k - 2} 𝐋^{(2)} 𝐑^{(2)} 𝐑^{(1)} = \dots = 𝐋^{(k)} \prod_{1}^{i = k} 𝐑^{(i)}

Věta (kritérium konvergence trojúhelníkové metody). Nechť je matice

𝐀

regulární a má všechna vlastní čísla jednonásobná a absolutně různá. Nechť existují LU rozklady matic

𝐗

𝐗^{- 1} 𝐋^{(0)}

{𝐀 𝐋}^{(k)}

od nějakého

k

. Potom posloupnosti

𝐋^{(k)}, 𝐑^{(k)}

konvergují a na diagonále

𝐑

je spektrum matice

𝐀

seřazené podle absolutní hodnoty.

Důkaz.

𝐀^{k} 𝐋^{(0)} = {𝐗 𝐃}^{k} 𝐗^{- 1} 𝐋^{(0)}

Nechť

𝐗 ≕ 𝐋_{𝐗} 𝐑_{𝐗}, 𝐗^{- 1} 𝐋^{(0)} ≔ 𝐋_{𝐘} 𝐑_{𝐘}

𝐀^{k} 𝐋^{(0)} = 𝐋_{𝐗} 𝐑_{𝐗} 𝐃^{k} 𝐋_{𝐘} 𝐑_{𝐘} = 𝐋_{𝐗} 𝐑_{𝐗} 𝐃^{k} 𝐋_{𝐘} 𝐃^{- k} 𝐃^{k} 𝐑_{𝐘}

{(𝐃^{k} {𝐋_{𝐘} 𝐃}^{- k})}_{i, j} = {\begin{matrix} 0, & i < j \\ 1, & i = j \\ λ_{i}^{k} L_{𝐘, i, j} λ_{j}^{- k}, & i > j \end{matrix}

Bez újmy na obecnosti můžeme předpokládat, že v

𝐃

jsou vlastní čísla seřázena absolutně sestupně, takže

| \frac{λ_{i}}{λ_{j}} | < 1

. Tedy nediagonální část této matice konverguje k nule a můžeme psát

𝐃^{k} {𝐋_{𝐘} 𝐃}^{- k} = 𝐈 + 𝐄^{(k)}, 𝐄^{(k)} \to 𝟎

\begin{aligned} 𝐀^{k} 𝐋^{(0)} & = 𝐋_{𝐗} 𝐑_{𝐗} (𝐈 + 𝐄^{(k)}) 𝐃^{k} 𝐑_{𝐘} \\ = 𝐋_{𝐗} 𝐑_{𝐗} (𝐈 + 𝐄^{(k)}) {𝐑_{𝐗}}^{- 1} 𝐑_{𝐗} 𝐃^{k} 𝐑_{𝐘} \\ = 𝐋_{𝐗} (𝐈 + 𝐑_{𝐗} 𝐄^{(k)} {𝐑_{𝐗}}^{- 1}) 𝐑_{𝐗} 𝐃^{k} 𝐑_{𝐘} \\ ≕ 𝐋_{𝐗} {𝐋_{𝐄}}^{(k)} {𝐑_{𝐄}}^{(k)} 𝐑_{𝐗} 𝐃^{k} 𝐑_{𝐘} \end{aligned}

LU rozklad můžeme provést, protože matice konverguje k identitě, z čehož zároveň víme, že také

{𝐋_{𝐄}}^{(k)}, {𝐑_{𝐄}}^{(k)} \to 𝐈

. Podařilo se nám tedy rozepsat

𝐀^{k} 𝐋^{(0)}

v LU tvaru a podle předchozího lemmatu víme, že levá část se rovná

𝐋^{(k)}

, tudíž

𝐋^{(k)} = 𝐋_{𝐗} {𝐋_{𝐄}}^{(k)} \to 𝐋_{𝐗}

. Nyní můžeme zlimitit iterační předpis a dostáváme

𝐑 = 𝐋^{- 1} 𝐀 𝐋 = {𝐋_{𝐗}}^{- 1} 𝐀 𝐋_{𝐱} = {𝐋_{𝐗}}^{- 1} {𝐗 𝐃 𝐗}^{- 1} 𝐋_{𝐱} = {𝐋_{𝐗}}^{- 1} {𝐋_{𝐗} 𝐑_{𝐗} 𝐃 𝐑_{𝐗}}^{- 1} {𝐋_{𝐗}}^{- 1} 𝐋_{𝐱} = {𝐑_{𝐗} 𝐃 𝐑_{𝐗}}^{- 1}

což je součin horních trojúhelníkových matic, takže diagonála se rovná součinu diagonál. Tento součin je zjevně roven diagonále

𝐃

, kde jsou vlastní čísla seřazená podle absolutní hodnoty, což mělo být dokázáno.

Pak je tu ještě LR algoritmus, který spočívá jednoduše v tom, že matici opakovaně rozložíme do LU rozkladu $𝐀^{(k)} ≕ {\hat{𝐋}}^{(k)} {\hat{𝐑}}^{(k)}$ a pak ho vynásobíme ve vopáčném pořadí: $𝐀^{(k + 1)} ≔ {\hat{𝐑}}^{(k)} {\hat{𝐋}}^{(k)}$ .

Věta. Všechny členy posloupnosti

𝐀^{(k)}

jsou si podobné.

Důkaz.

𝐀^{(k + 1)} = {\hat{𝐑}}^{(k)} {\hat{𝐋}}^{(k)} = {({\hat{𝐋}}^{(k)})}^{- 1} 𝐀^{(k)} {\hat{𝐋}}^{(k)} ■

LR algoritmus má menší nárok na paměť — stačí si pamatovat dvě matice a ne tři; nemusíme si pamatovat matici $𝐀$ . Ovšem tím, že si ji nepamatujeme, ztrácíme samoopravovací vlastnost metody. Algoritmus je obecně numericky nestabilní, proto se moc nepoužívá.

Věta. Existuje-li LU rozklad

𝐀^{k} ≕ ℒ^{(k)} ℛ^{(k)}

, potom

ℒ^{(k)} = \prod_{i = 1}^{k} {\hat{𝐋}}^{(i)}

ℛ^{(k)} = \prod_{1}^{i = k} {\hat{𝐑}}^{(i)}

Důsledek. Pokud v trojúhelníkové metodě zvolíme

𝐋^{(0)} ≔ 𝐈

, potom

𝐋^{(k)} = \prod_{i = 1}^{k} {\hat{𝐋}}^{(i)}

𝐑^{(k)} = {\hat{𝐑}}^{(k)}

Důsledek. Pokud konverguje trojúhelníková metoda začínající v jednotkové matici, potom LR algoritmus konverguje k horní trojúhelníkové matici.

Důkaz.

𝐀^{(k)} = {\hat{𝐋}}^{(k)} {\hat{𝐑}}^{(k)} = {(𝐋^{(k - 1)})}^{- 1} 𝐋^{(k)} 𝐑^{(k)} \to 𝐋^{- 1} 𝐋 𝐑 = 𝐑 ■

Obě tyto metody jsou náročné, protože každá iterace má složitost $Θ (n^{3})$ (kvůli výpočtu LU rozkladu). Také jsou poněkud numericky nestabilní, ale jsou dobrým teoretickým podkladem pro QR algoritmus, který už je lepší.

QR rozklad

Věta. Nechť

𝐀

je regulární matice. Potom existuje rozklad

𝐀 = 𝐐 𝐑

, kde

𝐐

je unitární a

𝐑

je horní trojúhelníková s kladnou diagonálou.

Důkaz (jednoznačnosti). Nechť

𝐀 = 𝐐_{1} 𝐐_{2}

. Z toho dostaneme:

𝐐_{2}^{T} 𝐐_{1} = 𝐑_{2} 𝐑_{1}^{- 1}

𝐐_{1}^{T} 𝐐_{2} = 𝐑_{1} 𝐑_{2}^{- 1}

Z toho plyne, že obě matice na levé straně jsou horní trojúhleníkové. Jelikož jsou vzájemně transponované, musí být diagonální, tedy

𝐐_{1}^{T} 𝐐_{2} ≕ 𝐃

. Dále máme:

𝐐_{2} = 𝐐_{1} 𝐐_{1}^{T} 𝐐_{2} = 𝐐_{1} 𝐃

𝐑_{2} = 𝐐_{2}^{T} 𝐐_{2} 𝐑_{2} = 𝐐_{2}^{T} 𝐐_{1} 𝐑_{1} = 𝐃 𝐑_{1}

𝐈 = 𝐐_{2}^{T} 𝐐_{2} = {(𝐐_{1} 𝐃)}^{T} 𝐐_{1} 𝐃 = 𝐃 𝐐_{1}^{T} 𝐐_{1} 𝐃 = 𝐃^{2}

Tedy diagonála

𝐃

obsahuje jen

\pm 1

. Jelikož

𝐑_{2} = 𝐃 𝐑_{1}

, a obě

𝐑

mají kladnou diagonálu, musí být

𝐃 = 𝐈

, tedy rozklady se rovnají.

QR rozklad pomocí Gram-Schmidtova ortogonalizačního procesu

Již známe z prváku:

\vec{{\tilde{q}}_{i}} ≔ {\vec{x}}_{i} - \sum_{j = 1}^{i - 1} ({\vec{x}}_{i} \cdot {\vec{q}}_{j}) {\vec{q}}_{j}

{\vec{q}}_{i} ≔ {\frac{\vec{{\tilde{q}}_{i}}}{‖ \vec{{\tilde{q}}_{i}} ‖}}_{2}

Když z těchto vzorců vyjádříme ${\vec{x}}_{i}$ pomocí ${\vec{q}}_{j}$ , dostaneme $𝐗 = 𝐐 𝐑$ , kde $𝐗$ se skládá z ortonormálníčh vektorů a $𝐑$ je horní trojůhelníková s kladnou diagonálou. Tedy Gram-Schmidt je algoritmus na výpočet QR rozkladu!

Ve skutečnosti se pro zvýšení stability používá mírná modifikace, kde při výpočtu ${\vec{q}}_{i}$ nepočítáme skalární součin s ${\vec{x}}_{i}$ , ale s již napočtenou verzí ${\vec{q}}_{i}$ , což samozřejmě v teorii neovlivní výsledek. I tak je ale dost nestabilní, takže se spíš používá v metodách pro řešení soustav rovnic. Každopádně má složitost $Θ (n^{3})$ , takže jsme si zatím moc nepomohli.

QR rozklad pomocí Householderových reflexí

Householderova reflexe nám dává způsob, jak pomocí unitární transformace zobrazit vektor na libovolný jiný vektor se stejnou normou. Takže speciálně pokud si zvolíme správný násobek jednotkového vektoru, můžeme ji použít k tomu, abychom všechny složky kromě jedné vynulovali. Akorát si musíme dávat pozor na numerickou stabilitu — nechceme, aby vektory byly příliš blízko, protože bychom potom dělili malým číslem, takže pokud příslušná složka vektoru bude kladná, budeme ho reflektovat do mínus jednotkového vektoru. Pokud chceme zachovat více složek, tak si Householderovu matici trochu upravíme, aby pár prvních složek ignorovala.

Tohle můžeme podobným způsobem jako u Gaussovy eliminační metody využít k převodu matice do trojůhelníkového tvaru. Nejdřív vynulujeme v prvním sloupci všechno kromě prvního prvku (ten bude nahrazen nějakým kladným číslem). Potom už první řádek necháme na pokoji a pomocí druhého řádku vynulujeme ve druhém sloupci další řádky. A tak dále.

Součinem všech Householderových matic, které jsme použili, dostaneme $𝐐$ . Ovšem kdybychom ten součin počítali jen tak, vyšla by z toho složitost $Θ (n^{4})$ . Naštěstí to můžeme spočítat efektivněji tak, že si je rozepíšeme z definice a budeme je na sebe navzájem aplikovat:

𝐇_{\vec{w}} 𝐀 = 𝐀 - 2 \vec{w} {\vec{w}}^{*} 𝐀 = 𝐀 - 2 \vec{w} {(𝐀^{*} \vec{w})}^{*}

QR rozklad pomocí Givensových rotací

Uděláme přesně to samé — budeme unitárně transformovat vektor do násobku jednotkoveho vektoru a tím nulovat složky — ale místo reflexe využijeme rotaci.

Definice. Givensova rotace je matice ve tvaru

G (i, j, θ) ≔ (\begin{array}{c} 1 \\ ⋱ \\ \cos (θ) & \dots & \sin (θ) \\ ⋮ & ⋱ & ⋮ \\ - \sin (θ) & \dots & \cos (θ) \\ ⋱ \\ 1 \end{array})

Jde tedy o unitární transformaci, která dokáže $j$ -tou složku vektoru nastavit na nulu, $i$ -tou na $\sqrt{x_{i}^{2} + x_{j}^{2}}$ a ty zbylé nechat. Takže přináší „jemnější“ způsob, jak dosáhnout toho samého. Inverzi Givensovy rotace získám jednoduše výměnou $\sin$ a $- \sin$ , tedy transpozicí.

Podobně jako předtím: Givensovy rotace nebudeme násobit jako idioti, ale využijeme toho, že se liší jen o pár složek od jednotkových matic. Složitost je potom opět $Θ (n^{3})$ .

QR algoritmus

Funguje úplně stejně jako LR algoritmus, akorát místo LU rozkladu počítáme QR rozklad. Posloupnosti se tentokrát budou jmenovat $𝐓^{(k)}$ , $𝐐^{(k)}$ a $𝐑^{(k)}$ .

Věta. Všechny členy posloupnosti

𝐓^{(k)}

jsou si ortogonálně podobné.

Důkaz.

𝐓^{(k + 1)} = 𝐑^{(k)} 𝐐^{(k)} = {(𝐐^{(k)})}^{- 1} 𝐀^{(k)} 𝐐^{(k)} ■

Pokud bude konvergovat $\prod_{i = 0}^{k} 𝐐^{(i)} \to 𝐔$ a $𝐓^{(k)} \to 𝐑$ , dostáváme Schurův rozklad: $𝐀 = 𝐔^{*} 𝐑 𝐔$ . Z diagonált $𝐑$ opět vyčteme spektrum $𝐀$ .

Lemma. Existuje-li QR rozklad

𝐀^{k} ≕ 𝒬^{(k)} ℛ^{(k)}

, potom

𝒬^{(k)} = \prod_{i = 1}^{k} {\hat{𝐐}}^{(i)}

ℛ^{(k)} = \prod_{1}^{i = k} {\hat{𝐑}}^{(i)}

Důkaz. Analogicky jako u LR algoritmu.

Věta. Pokud má matice

𝐀

absolutně různá vlastní čísla, potom posloupnost

𝐓^{(k)}

konverguje k horní trojúhelníkové matici s vlastními čísly seřazenými podle absolutní hodnoty na diagonále. Navíc je-li

𝐀

symetrická, potom limitní matice je diagonální.

Důkaz.

𝐀^{k} = {𝐗 𝐃}^{k} 𝐗^{- 1}

Nechť

𝐗 ≕ 𝐐_{𝐗} 𝐑_{𝐗}, 𝐗^{- 1} ≔ 𝐋_{𝐘} 𝐑_{𝐘}

𝐀^{k} = 𝐐_{𝐗} 𝐑_{𝐗} 𝐃^{k} 𝐋_{𝐘} 𝐑_{𝐘} = 𝐐_{𝐗} 𝐑_{𝐗} 𝐃^{k} 𝐋_{𝐘} 𝐃^{- k} 𝐃^{k} 𝐑_{𝐘}

Analogicky jako u LR algoritmu je

𝐃^{k} {𝐋_{𝐘} 𝐃}^{- k} = 𝐈 + 𝐅^{(k)}, 𝐅^{(k)} \to 𝟎

\begin{aligned} 𝐀^{k} & = 𝐐_{𝐗} 𝐑_{𝐗} (𝐈 + 𝐅^{(k)}) 𝐃^{k} 𝐑_{𝐘} \\ = 𝐐_{𝐗} 𝐑_{𝐗} (𝐈 + 𝐅^{(k)}) {𝐑_{𝐗}}^{- 1} 𝐑_{𝐗} 𝐃^{k} 𝐑_{𝐘} \\ = 𝐐_{𝐗} (𝐈 + 𝐑_{𝐗} 𝐅^{(k)} {𝐑_{𝐗}}^{- 1}) 𝐑_{𝐗} 𝐃^{k} 𝐑_{𝐘} \\ ≕ 𝐐_{𝐗} {𝐐_{𝐆}}^{(k)} {𝐑_{𝐆}}^{(k)} 𝐑_{𝐗} 𝐃^{k} 𝐑_{𝐘} \end{aligned}

QR rozklad můžeme provést, protože matice konverguje k identitě, z čehož zároveň víme, že také

{𝐐_{𝐆}}^{(k)}, {𝐑_{𝐆}}^{(k)} \to 𝐈

. Podařilo se nám tedy rozepsat

𝐀^{k}

v QR tvaru a podle předchozího lemmatu víme, že levá část se rovná

𝒬^{(k)} = \prod_{i = 1}^{k} 𝐐^{(i)}

, tudíž

𝐓^{(k)} = {𝒬^{(k)}}^{*} 𝐀 𝒬^{(k)} \to 𝐐_{𝐗}^{*} 𝐀 𝐐_{𝐗} ≕ 𝐓

. Máme

𝐓^{(k)} = {𝒬^{(k)}}^{*} {𝐗 𝐃 𝐗}^{- 1} 𝒬^{(k)} = {𝐐_{𝐆}}^{(k)}^{*} 𝐐_{𝐗}^{*} {𝐐_{𝐗} 𝐑_{𝐗} 𝐃 𝐑_{𝐗}}^{- 1} 𝐐_{𝐗}^{*} {𝐐_{𝐗} 𝐐_{𝐆}}^{(k)} \to {𝐑_{𝐗} 𝐃 𝐑_{𝐗}}^{- 1}

což je součin horních trojúhelníkových matic, takže diagonála se rovná součinu diagonál. Tento součin je zjevně roven diagonále

𝐃

, kde jsou vlastní čísla seřazená podle absolutní hodnoty. Navíc je-li

𝐀

symetrická, platí

𝐐_{𝐗}^{T} 𝐓 𝐐_{𝐗} = 𝐀 = 𝐀^{T} = 𝐐_{𝐗}^{T} 𝐓^{T} 𝐐_{𝐗} ∴ 𝐓 = 𝐓^{T}

a jelikož

𝐓

je trojúhelníková, musí být diagonální.

Zatím to vypadá, že QR algoritmus má stejně jako LR algoritmus složitost $Θ (n^{3})$ . Ale pomocí Hessenbergových QR iterací se dá zredukovat na $Θ (n^{2})$ .

Definice. Matice je v Hessenbergově tvaru, pokud je horní trojúhelníková s tím, že může mít také nenulové prvky na spodní diagonále.

Do Hessenbergova tvaru lze libovolnou matici převést podobnostní transformací přímým algoritmem se složitostí $Θ (n^{3})$ . Použijeme Housholderovy reflexe podobně jako pro výpočet QR rozkladu, ale s tím, že první řádek budeme už od začátku nechávat na pokoji. To nám zajistí, že když matici poté vynásobíme tou samou reflexí i zprava, abychom dostali podobnou matici, nerozbijeme si tím vytvořené nuly.

Vylepšený QR algoritmus bude fungovat tak, že $𝐀$ nejprve převedeme do Hessenbergova tvaru. To je sice $Θ (n^{3})$ , ale aspoň to budeme provádět jen jednou a ne při každé iteraci. Posloupnost $𝐓^{(k)}$ tentokrát budeme značit $𝐇^{(k)}$ . Nyní když chceme spočítat QR rozklad, tak nám stačí eliminovat prvky na spodní diagonále, a na to stačí $n - 1$ Givensových rotací! Ještě se potřebujeme ujistit, že po provedení zbytku kroku opět vznikne matice v Hessenbergově tvaru.

Věta. Součin

𝐑^{(k)} \prod_{i = 1}^{n - 1} 𝐆_{n - 1}^{(k)}

je opět matice v Hessenbergově tvaru.

Důkaz. Givensovy rotace, které používáme, vždy mění pouze dva po sobě jdoucí řádky, takže to takhle vyjde.

Nelineární rovnice

Začneme s funkcí jedné reálné proměnné. Nejprve musíme separovat kořeny — najít takové intervaly, aby v každém byl pouze jeden kořen. Metody totiž nemusí konvergovat, když máme kořenů víc. Na separaci není obecně žádná metoda, musí se to udělat pro konkrétní problém.

Nejjednodušší je bisekce — poněkud pomalá, ale velmi robustní. Odhad chyby můžeme provádět přímo pomocí residua $| f (x) |$ . Pokud ale funkce u kořene roste/klesá hodně pomalu, může tím vzniknout zbytečně špatná aproximace. Místo toho tedy můžeme použít odhad $| \frac{f (x)}{f^{'} (x)} |$ . Bisekce samozřejmě může naprosto selhat, pokud je funkce nespojitá.

Mějme nějaký odhad $x$ . Využitím Lagrangeovy věty o přírůstku dostáváme pro kořen:

α = x - \frac{f (x)}{f (ξ)}, ξ \in ⁄ α, x ⁄

Tohoto můžeme využít k vytvoření iterativní metody, ovšem musíme umět nějak odhadnout $f (ξ)$ .

Věta. Nechť

α \in ℝ, φ : ℝ \to ℝ

. Nechť platí

φ (α) = α

φ

je diferencovatelná na

V ≔ ⟨ α - r, α + r ⟩

\sup {| φ^{'} (x) | | x \in V} ≕ K < 1

. Potom posloupnost

x_{k} ≔ φ^{k} (x)

konverguje k

α

pro

x_{0} \in V

Poznámka. Pokud je

φ^{'}

spojitá v

α

, potom je podmínka splněna, pokud

φ^{'} (α) < 1

Důkaz. Použijeme větu o přírůstku funkce:

| φ (x_{k}) - α | = | φ (x_{k}) - φ (α) | = | φ^{'} (ξ_{k}) | | x_{k} - α | \leq K | x_{k} - α |

Tedy

| x_{k} - α | \leq K^{k} | x_{k} - α |

, z čehož plyne konvergence.

Definice. Iterativní metoda daná vztahem

x_{k} ≔ φ^{k} (x)

má řád konvergence

m

, pokud pro nějakou konstantu

C

platí

| x_{k + 1} - α | \leq C {| x_{k} - α |}^{m}

Věta. Nechť

φ \in 𝒞^{m} (H_{α})

(\forall i \in \hat{m - 1}) (φ^{(i)} (α) = 0)

. Potom

(\forall x_{k} \in H_{α}) (\exists ξ \in ⁄ α, x ⁄) (x_{k + 1} - α = \frac{φ^{(m)} (ξ)}{m!} {(x_{k} - α)}^{m})

Důkaz. Plyne přímo z Lagrangeova zbytku Taylorova rozvoje.

Metoda regula falsi funguje podobně jako bisekce, ale místo sekání napůl budeme sekat tam, kde se protíná úsečka mezi body s osou $x$ . Abychom to vyjádřili v obecném tvaru, máme

φ (x) ≔ \frac{x - x^{'}}{f (x) - f (x^{'})} f (x)

kde $x^{'}$ je druhý bod.

Věta. Nechť

f \in 𝒞^{1} (H_{α})

f^{'} (α) \neq 0

. Potom existuje takové

B_{α} \subseteq H_{α}

, že metoda regula falsi konverguje na

B_{α}

rychlostí prvního řádu.

Důkaz. Nechť

x_{k}, x_{k}^{'} \in B_{α}

. Použijeme dvakrát větu o přírůstku funkce:

\begin{aligned} | x_{k + 1} - α | & = | x_{k} - \frac{x_{k} - x_{k}^{'}}{f (x_{k}) - f (x_{k}^{'})} f (x_{k}) - α | \\ = | x_{k} - α - \frac{f (x_{k}) - f (α)}{f^{'} (ξ_{0})} | \\ = | x_{k} - α - \frac{f^{'} (ξ_{1}) (x_{k} - α)}{f^{'} (ξ_{0})} | \\ = | \frac{f^{'} (ξ_{1}) - f^{'} (ξ_{0})}{f^{'} (ξ_{1})} | | x_{k} - α | \end{aligned}

Jelikož

f^{'} (α) \neq 0

, pro

ξ

dostatečně blízké

α

existuje takové

C

, že

| f^{'} (ξ) | > C

, tedy

| \frac{f^{'} (ξ_{1}) - f^{'} (ξ_{0})}{f^{'} (ξ_{1})} | \leq \frac{1}{C} | f^{'} (ξ_{1}) - f^{'} (ξ_{0}) |

Jelikož derivace je spojitá, tento součin dokážeme udělat libovolně malý.

Pokud v metodě regula falsi pošleme $x^{'} \to x$ , dostáváme Newtonovu metodu:

φ (x) ≔ x - \frac{f (x)}{f^{'} (x)}

Věta. Nechť

f \in 𝒞^{1} (H_{α})

f^{'} (α) \neq 0

. Potom existuje takové

B_{α} \subseteq H_{α}

, že Newtonova metoda konverguje na

B_{α}

rychlostí prvního řádu.

Důkaz. Analogicky jako u metody regula falsi, akorát si ušetříme jedno použití věty o přírůstku funkce.

Věta. Nechť

f \in 𝒞^{2} (H_{α})

f^{'} (α) \neq 0

. Potom existuje takové

B_{α} \subseteq H_{α}

, že Newtonova metoda konverguje na

B_{α}

rychlostí druhého řádu.

Důkaz. Z důkazu předchozí véty máme

| x_{k + 1} - α | \leq \frac{1}{C} | f^{'} (x_{k}) - f^{'} (ξ_{0}) | | x_{k} - α |

Použijeme znovu větu o přírůstku funkce:

| x_{k + 1} - α | \leq \frac{1}{C} | f^{''} (ξ_{1}) | | x_{k} - ξ_{0} | | x_{k} - α | \leq \frac{1}{C} | f^{''} (ξ) | {| x_{k} - α |}^{2}

Obecně Newtonova metoda konverguje rychleji, ale vyžaduje přesnější odhad. Dá se udělat to, že pár kroků uděláme bisekcí nebo regula falsi a následně použijeme Newtonovu metodu.

Pomocí Čebyšenovy metody se dají odvodit podobné metody vyššího řádu, ale vyžadují ještě přesnější odhad a v praxi se moc nepoužívají.

Tyto metody často nekonvergují kvůli tomu, že dělají příliš velký krok. To se dá vylepšit tím, že se v každém kroku nejdřív podíváme, kam bychom skočili, a pokud by se tím residuum zvýšilo, tak zmenšíme skok (například na polovinu) a zkusíme znovu. Tomu se potom říká globálně konvergující metoda, protože pro každý počáteční odhad buď konverguje (i když hodně pomalu), nebo selže.

Nyní se podívejme na soustavy nelineárních rovnic. Ukazuje se, že stačí trochu zobecnit Newtonovu metodu:

φ (\vec{x}) ≔ \vec{x} - 𝐉_{\vec{f}}^{- 1} (\vec{x}) \vec{f} (\vec{x})

neboli v praktičtějším tvaru:

𝐉_{\vec{f}} ({\vec{x}}^{(k)}) ({\vec{x}}^{(k + 1)} - {\vec{x}}^{(k)}) = - \vec{f} ({\vec{x}}^{(k)})

Opět stačí při řešení lineární soustavy provést jen pár kroků iterativní metody, protože pro malá $k$ nám stačí malá přesnost a pro velká $k$ máme dobrý odhad.

Věta (o přírůstku funkce). Nechť

f \in 𝒞^{1} (H)

, kde

H

je konvexní oblast. Potom

(\forall \vec{u}, \vec{v} \in H) (\exists ξ \in H) (f (\vec{u}) - f (\vec{v}) = \vec{\nabla} f (\vec{ξ}) (\vec{u} - \vec{v}))

Důkaz. Viz skripta z ANA3. Vezmeme úsečku mezi body a aplikujeme na ni normální větu o přírůstku.

Věta. Nechť

\vec{f} \in 𝒞^{1} (H_{\vec{a}})

𝐉_{\vec{f}} (\vec{a})

je regulární. Potom existuje takové

B_{\vec{a}} \subseteq H_{\vec{a}}

, že zobecněná Newtonova metoda konverguje na

B_{a}

rychlostí prvního řádu.

Důkaz.

\begin{aligned} ‖ {\vec{x}}^{(k + 1)} - \vec{a} ‖ & = ‖ {\vec{x}}^{(k)} - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)}) (f ({\vec{x}}^{(k)}) - f (\vec{a})) - \vec{a} ‖ \\ = ‖ {\vec{x}}^{(k)} - \vec{a} - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)}) 𝐉_{\vec{f}} ({\vec{ξ}}_{1}, \dots, {\vec{ξ}}_{n}) ({\vec{x}}^{(k)} - \vec{a}) ‖ \\ \leq ‖ 𝐈 - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)}) 𝐉_{\vec{f}} ({\vec{ξ}}_{1}, \dots, {\vec{ξ}}_{n}) ‖ ‖ {\vec{x}}^{(k)} - \vec{a} ‖ \end{aligned}

kde značíme

𝐉_{\vec{f}} ({\vec{ξ}}_{1}, \dots, {\vec{ξ}}_{n}) ≔ (\begin{array}{c} \vec{\nabla} f_{1} (ξ_{1}) \\ ⋮ \\ \vec{\nabla} f_{n} (ξ_{n}) \end{array})

Chceme opět dokázat, že první normu můžeme omezit nějakou konstantou menší než

1

\begin{aligned} ‖ 𝐈 - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)}) 𝐉_{\vec{f}} ({\vec{ξ}}_{1}, \dots, {\vec{ξ}}_{n}) ‖ & = ‖ 𝐉_{\vec{f}}^{- 1} (\vec{a}) 𝐉_{\vec{f}} (\vec{a}) - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)}) 𝐉_{\vec{f}} ({\vec{ξ}}_{1}, \dots, {\vec{ξ}}_{n}) ‖ \\ = ‖ (𝐉_{\vec{f}}^{- 1} (\vec{a}) - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)})) 𝐉_{\vec{f}} (\vec{a}) - 𝐉_{\vec{f}}^{- 1} ({\vec{x}}^{(k)}) (𝐉_{\vec{f}} ({\vec{ξ}}_{1}, \dots, {\vec{ξ}}_{n}) - 𝐉_{\vec{f}} (\vec{a})) ‖ \end{aligned}

Oba rozdíly v zázvorkách dokážeme udělat libovolně malé, takže celé tuto dokážeme udělat libovolně nulové.

Lagrangeův polynom — konstrukce

Máme nějakou funkci a chceme ji aproximovat polynomem. Mohli bychom zkusit Taylora, ale ten potřebuje znát hodně derivací, což často nemáme k dispozici. Místo toho použijeme Lagrangeův polynom, který naopak vyžaduje hodnotu funkce ve více bodech. Dejme tomu, že známe funkční hodnoty v bodech $x_{0}, \dots, x_{n}$ a chceme na ně napasovat polynom $L_{n} (x) = \sum_{i = 0}^{n} a_{i} x^{i}$ . Poté prostě řešíme soustavu

(\begin{array}{ccc} x_{0}^{0} & \dots & x_{0}^{n} \\ ⋮ & ⋱ & ⋮ \\ x_{n}^{0} & \dots & x_{n}^{n} \end{array}) (\begin{array}{c} a_{0} \\ ⋮ \\ a_{n} \end{array}) = (\begin{array}{c} f (x_{0}) \\ ⋮ \\ f (x_{n}) \end{array})

𝐕_{x_{0}, \dots, x_{n}}^{(n)} \vec{a} = \vec{f}

kde $𝐕$ je Vandermondova matice.

Věta. Vandermondova matice je regulární.

Důkaz. Pro zjednodušení budu číslovat permutace od nuly.

\det (𝐕) = \sum_{π \in 𝕊_{n + 1}} sgn (π) \prod_{i = 0}^{n} V_{π (i), i} = \sum_{π \in 𝕊_{n + 1}} sgn (π) \prod_{i = 0}^{n} x_{π (i)}^{i}

Všimněme si, že jde o polynom více proměnných, kde každý člen je stupně

\frac{n (n + 1)}{2}

. Kdybychom všechny proměnné kromě

x_{i}

zafixovali, dostali bychom polynom jedné proměnné. Jelikož pro

x_{i} = x_{j}

je determinant nulový, tento polynom má kořeny

x_{j}

pro všechna

j \neq i

. Celkově tedy všechny výrazy ve tvaru

x_{i} - x_{j}, i \neq j

dělí determinant. Jelikož těchto výrazů je až na pořadí

\frac{n (n + 1)}{2}

, determinant se musí rovnat jejich součinu krát nějaká konstanta, tudíž žádné jiné kořeny mít nemůže.

Důsledek. Existuje právě jeden interpolační polynom.

Pokud budou body $x$ příliš blízko k sobě, Vandermondova matice může být špatně podmíněná. Naštěstí existují i jiné způsoby, jak Lagrangeův polynom spočítat.

Lagrangeův tvar: Definujeme polynomy

l_{i} ≔ \prod_{\begin{array}{c} j = 0 \\ j \neq i \end{array}}^{n} \frac{x - x_{j}}{x_{i} - x_{j}}

Zřejmě $l_{i} (x_{j}) = δ_{i j}$ . Pomocí nich můžeme vyjádřit interpolační polynom jako

l_{n} (x) ≔ \sum_{i = 0}^{n} f (x_{i}) l_{j} (x)

Tato konstrukce je však poněkud výpočetně náročná — při každé změně nějakého $x_{i}$ musíme celý polynom přepočítávat.

Lepší je v tomto Newtonova formule. Budeme postupovat indukcí. Víme, že existuje právě jeden polynom stupně $n - 1$ , který souhlasí v $n - 1$ bodech. Pomocí něj najdeme polynom stupně $n$ , který souhlasí ve všech $n$ bodech. Budeme ho hledat ve tvaru:

L_{n} (x) ≔ L_{n - 1} (x) + c_{n} \prod_{i = 0}^{n - 1} (x - x_{i})

Snadno si vyjádříme:

c_{n} ≔ \frac{f (x_{n}) - L_{n - 1} (x_{n})}{\prod_{i = 0}^{n - 1} (x_{n} - x_{i})}

Jednodušší způsob, jak to spočítat, jsou poměrné diference. Když si to rozepíšeme nerekurzivně, máme

L_{n} (x) = \sum_{i = 0}^{n} c_{i} \prod_{j = 0}^{i - 1} (x - x_{j})

Z toho dostaneme lineární soustavu

(\begin{array}{ccccc} 1 & 0 & 0 & \dots & 0 \\ 1 & (x_{1} - x_{0}) & 0 & \dots & 0 \\ 1 & (x_{2} - x_{0}) & (x_{2} - x_{0}) (x_{2} - x_{1}) & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & (x_{n} - x_{0}) & (x_{n} - x_{0}) (x_{n} - x_{1}) & \dots & \prod_{i = 0}^{n - 1} (x_{n} - x_{i}) \end{array}) (\begin{array}{c} c_{0} \\ ⋮ \\ c_{n} \end{array}) = (\begin{array}{c} f (x_{0}) \\ ⋮ \\ f (x_{n}) \end{array})

𝐁 \vec{c} = \vec{f}

Vidíme, že $c_{k}$ závisí jen na $f (x_{0}), \dots, f (x_{k})$ , což se z nějakého důvodu zapisuje jako $c_{k} = f [x_{0}, \dots, x_{k}]$ . Těmto hodnotám se říká poměrné diference.

Věta.

f [x_{i}, \dots, x_{i + k}] = \frac{f [x_{i + 1}, \dots, x_{i + k}] - f [x_{i}, \dots, x_{i + k - 1}]}{x_{i + k} - x_{i}}

Důkaz. Pro zjednodušení budeme počítat

f [x_{0}, \dots, x_{k}]

. Víme, že existuje právě jeden polynom

L ≔ L_{k} (x_{0}, \dots, x_{k}) (x) = f [x_{0}, \dots, x_{k}] \prod_{i = 0}^{k - 1} (x - x_{i}) + P_{k - 1} (x) = f [x_{0}, \dots, x_{k}] x^{k} + Q_{k - 1} (x)

a polynom

M ≔ L_{k - 1} (x_{0}, \dots, x_{k - 1}) (x) = f [x_{0}, \dots, x_{k - 1}] x^{k - 1} + Q_{k - 2} (x)

a také

N ≔ L_{k - 1} (x_{1}, \dots, x_{k}) (x) = f [x_{1}, \dots, x_{k}] x^{k - 1} + R_{k - 2} (x)

kde každý se v příslušných bodech rovná hodnotě

f

. Zkusíme vyjádřit ten první pomocí těch dalších dvou. Tedy máme dva polynomy, které se rovnají funkci v

k

bodech, a chceme, aby se rovnaly v

k + 1

bodech. Víme, že v bodech

x_{1}, \dots, x_{k - 1}

se polynomy

M, N

rovnají. Použijeme nějakou jejich afinní kombinaci:

(α M + (1 - α) N) (x_{i}) = f (x_{i})

Tato kombinace bude ještě navíc záviset na

x

. V bodě

x_{0}

se zjevně hodí mít

α = 1

a v bodě

x_{k}

α = 0

. Abychom dostali něco hezkého, definujeme tedy

α (x) ≔ \frac{x - x_{0}}{x_{k} - x_{0}}

Jelikož to je polynom stupně

1

, máme nový polynom stupně

k

\begin{aligned} S (x) & ≔ α (x) M (x) + (1 - α (x)) N (x) \\ = N (x) + α (x) (M (x) - N (x)) \\ = N (x) + \frac{x - x_{0}}{x_{k} - x_{0}} (M (x) - N (x)) \\ = N (x) + \frac{x - x_{0}}{x_{k} - x_{0}} (f [x_{0}, \dots, x_{k - 1}] x^{k - 1} + Q_{k - 2} (x) - f [x_{1}, \dots, x_{k}] x^{k - 1} - R_{k - 2} (x)) \\ = \frac{f [x_{i + 1}, \dots, x_{i + k}] - f [x_{i}, \dots, x_{i + k - 1}]}{x_{i + k} - x_{i}} x^{k} + T_{k - 2} (x) \end{aligned}

To se podle jednoznačnosti musí po koeficientech rovnat

L

, čímž je důkaz hotov.

Takhle si tedy můžeme napočítat postupně všechny poměrné diference se složitostí $Θ (n^{2})$ , přičemž poté využijeme jenom ty, které začínají na $x_{0}$ .

Lagrangeův polynom — analýza

Takže už umíme efektivně najít Lagrangeův polynom, ale jak zjistíme, jestli je dobrou aproximací?

Věta. Nechť

I_{x} \subseteq D_{f}

je nejmenší interval takový, že

{x, x_{0}, \dots, x_{n}} \subseteq I_{x}

f

má na něm derivaci řádu

n + 1

. Nechť

L_{n}

je její Lagrangeův polynom. Definujme

ω_{n} (x) ≔ \prod_{i = 0}^{n} (x - x_{i})

Potom existuje takové

ξ \in I_{x}

, že

R_{n} (x) ≔ f (x) - L_{n} (x) = \frac{f^{(n + 1)} (ξ)}{(n + 1)!} ω_{n} (x)

Důkaz. Definujme pomocnou funkci

Q (t) ≔ ω_{n} (x) R_{n} (t) + ω_{n} (t) R_{n} (x)

Tedy má

n + 2

kořenů. Zřejmě

Q (x_{i}) = Q (x) = 0

. Podle Rolleovy věty má funkce na každém intervalu mezi body nulovou derivaci, takže její derivace má

n - 1

kořenů. Pokud ji budeme derivovat dál, dostaneme se k tomu, že derivace řádu

n + 1

má jeden kořen

ξ

na intervalu

I_{x}

. Nyní si tuto derivaci spočteme:

R_{n}^{(n + 1)} = f^{(n + 1)} - L_{n}^{(n + 1)} = f^{(n + 1)}

ω_{n}^{(n + 1)} = (n + 1)!

Q^{(n + 1)} (t) = ω_{n} (x) R_{n}^{(n + 1)} (t) - ω_{n}^{(n + 1)} (t) R_{n} (x) = ω_{n} (x) f^{(n + 1)} (t) - (n + 1)! R_{n} (x)

Pokud dosadíme

t ≔ ξ

, dostáváme tvrzení věty.

Všimněme si, že $ω_{n}$ mimo interval $(x_{0}, x_{n})$ roste hodně rychle, takže Lagrangeův polynom se nehodí pro extrapolaci, pouze pro interpolaci.

Také vidíme, že k výpočtu nepotřebujeme znát žádné derivace, ale aby šlo o dobrou aproximaci, musí $n + 1$ . derivace existovat a být nízká.

Věta.

R_{n} (x) = f [x_{0}, \dots, x_{n}, x] ω_{n} (x)

Důkaz. Stačí použít Newtonovu formuli s tím, že chceme jakoby zkonstruovat polynom

L_{n + 1}

x_{n + 1} ≔ x

Důsledek.

f [x_{0}, \dots, x_{n}, x] = \frac{f^{(n + 1)} (ξ)}{(n + 1)!}

Důsledek. Lagrangeův polynom se dá zapsat ve tvaru

L_{n} (x) = \sum_{i = 0}^{n} \frac{f^{(i)} (ξ_{i})}{i!} \prod_{j = 0}^{i - 1} (x - x_{j})

což je velmi podobné Taylorovu rozvoji.

Definice. Mějme

x_{0} \in ℝ, H_{x_{0}}

a funkce

f, g : H_{x_{0}} \to ℝ

. Řekneme, že

f

aproximuje

g

H_{x_{0}}

s přesností řádu

r

, pokud

\lim_{x \to x_{0}} \frac{| f (x) - g (x) |}{{| x - x_{0} |}^{r}} = C > 0

Pro Lagrangeův polynom máme

\lim_{x \to x_{0}} \frac{| f (x) - L_{n} (x) |}{| x - x_{0} |} = \lim_{x \to x_{0}} \frac{| \frac{f^{(n + 1)} (ξ)}{(n + 1)!} ω_{n} (x) |}{| x - x_{0} |} = \lim_{x \to x_{0}} | \frac{\frac{f^{(n + 1)} (ξ)}{(n + 1)!}}{x - x_{0}} \prod_{i = 0}^{n} (x - x_{i}) | = | \frac{f^{(n + 1)} (ξ)}{(n + 1)!} \prod_{i = 1}^{n} (x_{0} - x_{i}) |

Tedy Lagrangeův polynom na okolí bodů $x_{i}$ aproximuje funkci s přesností prvního řádu. Toto nezávisí na volbě $n$ a obecně neplatí, že by nám více bodů dávalo lepší aproximaci! Pomůže nám to, pokud všechny vyšší derivace jsou blízké nule, ale jinak dochází k Rungovu jevu.

Tento problém řeší interpolace po částech. Rozkrájíme interval, na kterém chceme funkci aproximovat, a v každém kousku si vytvoříme jiný Lagrangeův polynom. Speciálně pokud na každém intervalu použijeme pouze jeho krajní body, dostaneme aproximaci lomenou čarou. Nevýhoda je, že nemůžeme zajistit, aby výsledná aproximace byla diferencovatelná. To řeší Hermitův polynom, což je v určitém smyslu kombinace Taylora a Lagrange, ale opět k němu potřebujeme znát derivace $f$ .

Lagrangeův polynom se dá zobecnit do více dimenzí. Nejsnadněji se zkonstruuje, pokud známe hodnoty bodů na nějaké mřížce, například v $ℝ^{2}$ budeme znát $f (x_{1}, y_{1}), f (x_{1}, y_{2}), f (x_{2}, y_{1}), f (x_{2}, y_{2})$ .

Aproximace derivace

Hodí se pro výpočet Newtonovy metody a pro řešení diferenciálních rovnic.

Máme diferencovatelnou funkci $f$ , kterou známe jen v konečném počtu bodů $x_{0}, \dots, x_{n}$ , a chceme zjisit její derivaci.

Již víme, že $f (x) = L_{n} (x) + R_{n} (x)$ . Pokud všechny tyhle věci nějak vyjádříme, můžeme vztah zderivovat. Zderivujme si bazické polynomy:

l_{j}^{'} (x) = \sum_{i \neq j} \frac{1}{x - x_{i}} \prod_{k \neq j} \frac{x - x_{k}}{x_{j} - x_{k}}

L_{n}^{'} (x) = \sum_{j = 0}^{n} f (x_{j}) l_{j}^{'} (x)

Ale prakticky je lepší derivovat konkrétní polynom.

R_{n} (x) = \frac{f^{(n + 1)} (ξ (x))}{(n + 1)!} ω_{n} (x)

Bohužel nevíme, jak přesně závisí $ξ$ na $x$ . Budeme předpokládat, že je aspoň diferencovatelné.

R_{n}^{(k)} (x) = \frac{\sum_{i = 0}^{k} (\binom{k}{i}) {(f^{(n + 1)} (ξ (x)))}^{(k - i)} ω_{n}^{(i)} (x)}{(n + 1)!}

Z toho vidíme, že $f$ musí mít derivaci řádu $n + k + 1$ , což není úplně příjemné. Zkusme si obecně zderivovat $ω$ :

ω_{n}^{(k)} (x) = \sum_{i_{1} = 0}^{n} \sum_{\begin{array}{c} i_{2} = 0 \\ i_{2} \neq i_{1} \end{array}}^{n} \dots \sum_{\begin{array}{c} i_{k} = 0 \\ i_{k} \neq i_{1} \\ ⋮ \\ i_{k} \neq i_{k - 1} \end{array}}^{n} \prod_{\begin{array}{c} j = 0 \\ j \neq i_{1} \\ ⋮ \\ j \neq i_{k} \end{array}}^{n} (x - x_{j})

Z tohohle výrazu je naprosto evidentní, že derivace $ω$ nemusí být obecně nulová, takže výsledek nebude přesný. Což netuším, proč by někdo očekával. Co kdybychom zkusili přibližovat body k sobě, abychom napodobili limitu v definici derivace? Zvolíme nějaké $h \in ℝ^{+}$ a budeme brát bydy $x_{i} ≔ x_{0} + i h, h \in {- m_{1}, \dots, m_{2}}$ . Nechť $n ≔ m_{1} + m_{2}, f_{i} ≔ f (x_{i})$ . Zkonstruujeme Lagrangeův polynom. Pokud vezmeme $x ≔ x_{0} + t h$ , můžeme ho vyjádřit jako

L_{n} (t) ≔ \sum_{i = - m_{1}}^{m_{2}} f_{i} l_{i} (t)

l_{i} (t) ≔ \prod_{\begin{array}{c} j = - m_{1} \\ j \neq i \end{array}}^{m_{2}} \frac{t - j}{t - i}

Pokud vyjádříme předchozí výraz pro $ω$ v závislosti na $t$ , dostaneme

ω_{n}^{(k)} (t) = h^{n + 1 - i} \sum_{i_{1} = - m_{1}}^{m_{2}} \sum_{\begin{array}{c} i_{2} = - m_{1} \\ i_{2} \neq i_{1} \end{array}}^{m_{2}} \dots \sum_{\begin{array}{c} i_{k} = - m_{1} \\ i_{k} \neq i_{1} \\ ⋮ \\ i_{k} \neq i_{k - 1} \end{array}}^{m_{2}} \prod_{\begin{array}{c} j = - m_{1} \\ j \neq i_{1} \\ ⋮ \\ j \neq i_{k} \end{array}}^{m_{2}} (t - j)

To pro $h \to 0$ jde k nule! Tedy čím menší zvolíme $h$ , tím lepší máme aproximaci, konkrétně s řádem $n + 1 - k$ .

Pro aproximaci derivace v některém z uzlů $x_{0}, \dots, x_{n}$ se používají konečné diference.

Věta. Nechť

f \in 𝒞^{3} ⟨ x_{0}, x_{1} ⟩, h ≔ x_{1} - x_{0}

. Potom

| \frac{f (x_{1}) - f (x_{0})}{h} - f^{'} (x_{0}) | = 𝒪 (h)

tedy dopředná konečná diference aproximuje derivaci s přesností prvního řádu.

Důkaz.

\begin{aligned} f (x) & = L_{1} (x) + R_{1} (x) \\ = f (x_{0}) + \frac{f (x_{1}) - f (x_{0})}{x_{1} - x_{0}} (x - x_{0}) + \frac{f^{''} (ξ (x))}{2} (x - x_{0}) (x - x_{1}) \end{aligned}

\begin{aligned} f^{'} (x) & = L_{1}^{'} (x) + R_{1}^{'} (x) \\ = \frac{f (x_{1}) - f (x_{0})}{x_{1} - x_{0}} + \frac{f^{'''} (ξ (x)) ξ^{'} (x)}{2} (x - x_{0}) (x - x_{1}) + \frac{f^{''} (ξ (x))}{2} ((x - x_{0}) + (x - x_{1})) \end{aligned}

f^{'} (x_{0}) = \frac{f (x_{1}) - f (x_{0})}{h} + \frac{f^{''} (ξ (x_{0}))}{2} h

f^{'} (x_{0}) - \frac{f (x_{1}) - f (x_{0})}{h} = \frac{f^{''} (ξ (x_{0}))}{2} h = 𝒪 (h)

Poznámka. Analogicky můžeme udělat zpětnou konečnou diferenci s uzly

x_{- 1}, x_{0}

Věta. Nechť

f \in 𝒞^{4} ⟨ x_{- 1}, x_{1} ⟩, h ≔ x_{1} - x_{0} = x_{0} - x_{- 1}

. Potom

| \frac{f (x_{1}) - f (x_{- 1})}{2 h} - f^{'} (x_{0}) | = 𝒪 (h^{2})

tedy centrální konečná diference aproximuje derivaci s přesností druhého řádu.

Důkaz.

\begin{aligned} f (x) & = L_{2} (x) + R_{2} (x) \\ = f (x_{- 1}) + f [x_{- 1}, x_{0}] (x - x_{- 1}) + f [x_{- 1}, x_{0}, x_{1}] (x - x_{- 1}) (x - x_{0}) + \frac{f^{'''} (ξ (x))}{6} (x - x_{- 1}) (x - x_{0}) (x - x_{1}) \end{aligned}

\begin{aligned} f^{'} (x) & = L_{1}^{'} (x) + R_{1}^{'} (x) \\ = f [x_{- 1}, x_{0}] + f [x_{- 1}, x_{0}, x_{1}] ((x - x_{- 1}) + (x - x_{0})) + \frac{f^{''''} (ξ (x)) ξ^{'} (x)}{6} (x - x_{- 1}) (x - x_{0}) (x - x_{1}) \\ + \frac{f^{'''} (ξ (x))}{6} ((x - x_{0}) (x - x_{1}) + (x - x_{- 1}) (x - x_{1}) + (x - x_{- 1}) (x - x_{0})) \end{aligned}

f^{'} (x_{0}) = \dots = \frac{f (x_{1}) - f (x_{- 1})}{2 h} + \frac{f^{'''} (ξ (x_{0}))}{6} h^{2}

Vidíme, že pomocí Lagrange se to dokazuje pěkně blbě a navíc je potřeba dost vysoká trída diferencovatelnosti. Co takhle použít Taylora?

Věta. Nechť

f \in 𝒞^{2} ⟨ x_{0}, x_{1} ⟩, h ≔ x_{1} - x_{0}

. Potom

| \frac{f (x_{1}) - f (x_{0})}{h} - f^{'} (x_{0}) | = 𝒪 (h)

tedy dopředná konečná diference aproximuje derivaci s přesností prvního řádu.

Důkaz.

f (x_{1}) = f (x_{0}) + f^{'} (x_{0}) (x_{1} - x_{0}) + \frac{f^{''} (ξ)}{2} {(x_{1} - x_{0})}^{2} = f (x_{0}) + h f^{'} (x_{0}) + \frac{h^{2}}{2} f^{''} (ξ)

\frac{f (x_{1}) - f (x_{0})}{h} - f^{'} (x_{0}) = \frac{f^{''} (ξ)}{2} h = 𝒪 (h)

Poznámka. Analogicky můžeme udělat zpětnou konečnou diferenci s uzly

x_{- 1}, x_{0}

Věta. Nechť

f \in 𝒞^{3} ⟨ x_{- 1}, x_{1} ⟩, h ≔ x_{1} - x_{0} = x_{0} - x_{- 1}

. Potom

| \frac{f (x_{1}) - f (x_{- 1})}{2 h} - f^{'} (x_{0}) | = 𝒪 (h^{2})

tedy centrální konečná diference aproximuje derivaci s přesností druhého řádu.

Důkaz.

f (x_{1}) = f (x_{0}) + f^{'} (x_{1}) (x_{1} - x_{0}) + \frac{f^{''} (x_{1})}{2} {(x_{1} - x_{0})}^{2} + \frac{f^{'''} (ξ_{1})}{6} {(x_{1} - x_{0})}^{3}

f (x_{- 1}) = f (x_{0}) + f^{'} (x_{- 1}) (x_{- 1} - x_{0}) + \frac{f^{''} (x_{- 1})}{2} {(x_{- 1} - x_{0})}^{2} + \frac{f^{'''} (ξ_{- 1})}{6} {(x_{- 1} - x_{0})}^{3}

f (x_{1}) - f (x_{- 1}) = 2 h f^{'} (x_{0}) + \frac{h^{3}}{6} (f^{'''} (ξ_{1}) + f^{'''} (ξ_{- 1}))

a tak dále.

Věta. Nechť

f \in 𝒞^{3} ⟨ x_{- 1}, x_{1} ⟩, h ≔ x_{1} - x_{0} = x_{0} - x_{- 1}

. Potom

| \frac{f (x_{1}) - 2 f (x_{0}) + f (x_{- 1})}{h^{2}} - f^{''} (x_{0}) | = 𝒪 (h)

tedy jiná centrální konečná diference aproximuje druhou derivaci s přesností prvního řádu.

Důkaz. Opět vyjádříme

f (x_{1}), f (x_{- 1})

pomocí Taylora druhého stupně, ale tentokrát je sečteme.

Věta. Nechť

f \in 𝒞^{4} ⟨ x_{- 1}, x_{1} ⟩, h ≔ x_{1} - x_{0} = x_{0} - x_{- 1}

. Potom

| \frac{f (x_{1}) - 2 f (x_{0}) + f (x_{- 1})}{h^{2}} - f^{''} (x_{0}) | = 𝒪 (h^{2})

tedy jiná centrální konečná diference aproximuje druhou derivaci s přesností druhého řádu.

Důkaz. V důkazu předchozí věty dostáváme

| \frac{f (x_{1}) - 2 f (x_{0}) + f (x_{- 1})}{h^{2}} - f^{''} (x_{0}) | = \frac{h}{6} | f^{'''} (ξ_{1}) - f^{'''} (ξ_{- 1}) |

Použijeme větu o přírůstku funkce:

f^{'''} (ξ_{1}) - f^{'''} (ξ_{- 1}) = f^{''''} (ξ) (ξ_{1} - ξ_{- 1})

Jelikož

ξ_{1} - ξ_{- 1} < 2 h

, můžeme z odhadu vysáknout další

h

Důkaz (jiný). Použijeme Taylorovy rozvoje třetího stupně a sečteme je.

Věta. Nechť

f \in 𝒞^{n + 1} (x_{0}, x_{n})

. Potom pro libovolné

k \in \hat{n}

existuje konečná diference pro aproximaci

k

-té derivace s přesností řádu

n + 1 - k

Důkaz. Rozepíšeme si Taylorův rozvoj do řádu

n

f (x) = \sum_{j = 0}^{n} \frac{f^{(j)} (x_{0})}{j!} {(x - x_{0})}^{j} + \frac{f^{(n + 1)} (ξ (x))}{(n + 1)!} {(x - x_{0})}^{n + 1}

Když do něj postupně dosadíme body

x_{- m_{1}}, \dots, x_{- 1}, x_{1}, \dots, x_{m_{2}}

, dostáváme

n

rovnic tvaru

f_{i} = f_{0} + \sum_{j = 1}^{n} C_{i, j} h^{j} f^{(j)} (x_{0}) + d_{i} h^{n + 1} f^{(n + 1)} (ξ_{i})

Tyto rovnice lineárně nakombinujeme tak, aby

k

-tá derivace zůstala právě jednou a ostatní kromě nulté a

n - 1

vypadly:

\sum_{i = - m_{1}}^{m_{2}} α_{i} f_{i} = \sum_{i = - m_{1}}^{m_{2}} α_{i} f_{0} + h^{k} f^{(k)} (x_{0}) + \sum_{i = - m_{1}}^{m_{2}} α_{i} d_{i} h^{n + 1} f^{(n + 1)} (ξ_{i})

To potom vydělíme

h_{k}

a máme návod pro spočtení

k

-té derivace s chybou

𝒪 (h^{n + k - 1})

Výpisky z Numerické matematiky 1

Kontakt

Zdroje

Organizace

Cíl

Matice

Rozklady matic

Posloupnosti, normy

Numerické výpočty

LU rozklad

Modifikace Gaussovy eliminační metody

Stacionární iterativní metody

Vlastní čísla

Mocninná metoda

Trojúhelníková metoda, LR algoritmus

QR rozklad

QR rozklad pomocí Gram-Schmidtova ortogonalizačního procesu

QR rozklad pomocí Householderových reflexí

QR rozklad pomocí Givensových rotací

QR algoritmus

Nelineární rovnice

Lagrangeův polynom — konstrukce

Lagrangeův polynom — analýza

Aproximace derivace

Numerický výpočet integrálu

Otázky ke zkoušce