2. Numerické metody (90%)

tags: matika, MA018

• Iterativní metody pro řešení nelineárních rovnic
  • Newtonova metoda a její modifikace.
• Přímé metody pro řešení systému lineárních rovnic
  • Gaussova eliminace
  • Jacobi
  • Gauss-Seidel
  • relaxační metody
• Numerická diferenciace, diferenciační schémata.

Keď tak, tu sú moje (Megi) poznámky, čo som si na skúšku písala
https://www.dropbox.com/s/mw22d9gqlljd8lo/Numeriky.pdf?dl=0

Chyby

• \(x\): přesná hodnota
• \(\tilde{x}\): aproximace \(x\)
• \(x - \tilde{x}\): absolutní chyba \(\tilde{x}\)
• \(|x - \tilde{x}| \le \alpha\): odhad absolutní chyby
• \(\displaystyle \frac{x - \tilde{x}}{x}\): relativní chyba
• \(\displaystyle |\frac{x - \tilde{x}}{x}| \le \delta\): odhad relativní chyby

Parametrická spojitost funkce
Zápis \(C^2[a,b]\) značí že interval \([a,b]\) je \(C^2\) spojitý.
Neplést se geometickou spojitostí \(G\), která se definuje pomocí tečných vektorů.

• \(C^{-}\) – Dva segmenty nenavazují
• \(C^{0}\) – Dva segmenty na sebe navazují
• \(C^{1}\) – Navazují na sebe první derivace funkce
• \(C^{2}\) – Navazují na sebe druhé derivace funkce.

Iterativní metody pro řešení nelineárních rovnic

Prostá iterační metoda (metoda pevného bodu)

• Používá se pro rovnice typu \(x = g(x)\)
• Funkce G je spojitá na intervalu \(I = [a, b]\)
• Řešení rovnice \(\xi\) (ksí) = pevný bod funkce \(g\)

Proces

• Zvolíme iterační funkci \(g\)
• Zvolíme počáteční hodnotu \(x_0\) náležící do intervalu \(I\)
• Položíme \(x_1 = g(x_0)\) obecně \(x_{k+1} = g(x_k)\)

Jestliže spojitá funkce \(g\) zobrazuje interval \(I\) do sebe,
tj. pro každé \(x \in I\) platí \(g(x) \in I\), pak na intervalu \(I\) existuje
alespoň jeden pevný bod \(\xi\) funkce \(g\).

Klasifikace pevných bodů

Pevný bod \(\xi\) funkce \(g\) se nazývá

• přitahující (atraktivní) pevný bod, jestliže existuje takové okolí \(V\) tohoto bodu \(\xi\), že pro každou počáteční aproximaci \(x_0 \in V\) posloupnost iterací \(\{x_k\}^\infty_{k = 0}\) konverguje k bodu \(\xi\).
• odpuzující (repulzivní) pevný bod, jestliže existuje takové okolí \(U\) bodu \(\xi\), že pro každou počáteční aproximaci \(x_0 \in U\), \(x_0 \neq \xi\), existuje takové \(k\), že \(x_k \notin U\).

Určení

Necht \(g \in C[a,b]\), \(g:[a; b]→[a,b]\) a nechť \(g\) má v bodě \(\xi\) derivaci.

• Je-li \(|g'(\xi)| < 1\), pak \(\xi\) je pritahující pevný bod.
• Je-li \(|g'(\xi)| > 1\), pak \(\xi\) je odpuzující pevný bod.

Newtonova metoda (Metoda tečen)

Metoda druhého řádu pro jednoduchý kořen \(\xi\)

Iterační funkce:
\begin{align*} g(x_{k+1}) = x_k - \frac{f(x_k)}{f'(x_k)} \end{align*}

Proces

• Uvažujeme rovnici \(f(x) = 0\)
• Zvolíme počáteční hodnotu \(x_0\)
• \(x_1\) = průsečík tečny k \(f\) v bodě \(x_0\) a osy \(x\)
• obecně: \(x_{k+1}\) je průsečík tečny k funkci \(f\) v bodě \(x_k\) s osou \(x\)

Nechť \(f \in C^2[a, b]\) viz Spojitost funkce . Nechť \(\xi \in [a,b]\) je kořenem rovnice \(f(x) = 0\) a \(f'(\xi) \neq 0\). Pak existuje \(\delta > 0\) tak, že posloupnost \(\{x_k\}^\infty_{k = 0}\) generovaná Newtonovou metodou konverguje k bodu \(\xi\) pro každou pocátecní aproximaci \(x_0 \in [\xi - \delta, \xi + \delta] \subseteq [a,b]\).

Fourierovy podmínky

• Nechť \(f \in C^2[a, b]\) a rovnice \(f(x) = 0\) má v intervalu jediný koren \(\xi\)
• Nechť \(f"\) a \(f'\) nemění znaménka na intervalu \([a,b]\) pričemž \(f'(x) \neq 0\) a \(∀ \in [a,b]\)

Nechť počáteční aproximace \(x_0\) je ten z krajních bodu \([a,b]\), v nemž znaménko funkce je stejné jako znaménko \(f''(x)\) na intervalu \([a,b]\). Pak posloupnost \(\{x_k\}^\infty_{k = 0}\) určená Newtonovou metodou konverguje monotonně k bodu \(\xi\).

Problémy

• \(\xi\) – inflexní bod, např. funkce \(arctan(x)\)
• Nevhodná volba iterace: \(f'(x_k) = 0\).
• Vícenásobný koren: \(f(\xi) = f'(\xi)\) → pomalá konvergence

Metoda sečen

metoda je řádu \((1+\sqrt{5})/2 \doteq 1,618\).

Iterační funkce:
\begin{align*} x_{k+1} = x_k - \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})} f(x_k) \end{align*}

Metoda regula falsi

metoda je řádu 1.

Iterační funkce:
\begin{align*} x_{k+1} = x_k - \frac{x_k - x_s}{f(x_k) - f(x_s)} f(x_k) \end{align*}
\(s = s(k)\) je nejvetší index takový, že \(f(x_k)f(x_s) < 0\), přitom \(f(x_0)f(x_1) < 0\)

Přímé metody pro řešení systému lineárních rovnic

Gaussova eliminace

Rovnice se přepíše do maticového zápisu
Jsou povoleny následující operace, které nezmění hodnost soustavy:

• násobení či dělení jednotlivých řádků nenulovým číslem
• prohazování libovolných řádků soustavy
• přičítání násobků jednotlivých řádků k jinému řádku

Příklad:
\begin{align*} \begin{array}{cc} 8x & -y & -2z & = 0 \\ -x & +7y & -z & = 10 \\ -2x & -y & +9z & = 23 \\ \end{array} \end{align*}
přepíšeme na:
\begin{align*} \left(\begin{array}{cc} 8 & -1 & -2 & 0 \\ -1 & 7 & -1 & 10 \\ -2 & -1 & 9 & 23 \\ \end{array}\right) \end{align*}

• První krok – eliminace x z druhého a třetího řádku.
  1. první řádek neupravujeme
  2. druhý řádek násobíme 8 a přičteme první řádek
  3. třetí řádek násobíme 4 a přičteme první řádek
     \begin{align*} \left(\begin{array}{cc} 8 & -1 & -2 & 0 \\ 0 & 55 & -10 & 80 \\ 0 & -5 & 34 & 92 \\ \end{array}\right) \end{align*}
• Druhý krok – eliminace y ze třetího řádku
  1. první řádek neupravujeme
  2. druhý řádek neupravujeme
  3. třetí řádek násobíme 11 a přičteme druhý řádek
  4. vydělíme třetí řádek 364 a získáme řešení z = 3
     \begin{align*} \left(\begin{array}{cc} 8 & -1 & -2 & 0 \\ 0 & 55 & -10 & 80 \\ 0 & 0 & 1 & 3 \\ \end{array}\right) \end{align*}
• Třetí krok – zpětné dosazení
  1. k druhému řádku přičteme třetí řádek vynásobený 10
  2. druhý řádek podělíme 55
  3. k prvímu řádku přičteme třetí řádek vynásobený 2 a druhý řádek

\begin{align*} \left(\begin{array}{cc} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 2 \\ 0 & 0 & 1 & 3 \\ \end{array}\right) \end{align*}
Výsledek je pak \(x = 1\), \(y = 2\), \(z = 3\)

Jacobiova iterační metoda

matici A zapíšeme jako
\(A = D + L + U\)

• D je diagonální matice
• L je dolní trojúhelníková matice s nulami na diagonále
• U je horní trojúhelníková matice s nulami na diagonále

\(Ax = (D + L + U)x = b\)
\(Dx = -(L + U)x + b\)

Pokud \(a_{ii} \neq 0, i = 1,...,n\), je matice \(D\) regulární a z předchozí rovnice lze vypočítat

\(x = \frac{-D(L + U)x + b}{D}\)

Jacobiova iterační matice: \(T_j = -D^{-1}(L + U)\)
\(\displaystyle T_J = (t_{ij}), t_{ij} = -\frac{a_{ij}}{a_{ii}}\) pro \(i \neq j, t_{ii} = 0\) pro \(i = 1,...,n\)

Výsledný maticový zápis:

Věta o konvergenci

Posloupnost \(\{x_k\}^\infty_{k = 0}\) generovaná metodou \(x^{k+1} = T_Jx^k + D^{-1}b\) konverguje pro každou počáteční aproximaci \(x^0 \in \mathbb{R}^n\) právě tehdy, když \(\varrho(T_J) < 1\).

Silné řádkové sumační kriterium

Pokud je matice ryze řádkově diagonálně dominantní pak Jacobiova iterační metoda konverguje pro každou počáteční aproximaci \(x^0 \in \mathbb{R}^n\)

\begin{align*} |a_{ii}| > \sum_{j = 1, j \neq i}^{n} |a_{ij}| ,\quad i = 1,...,n \end{align*}

Silné sloupcově sumační kriterium

Pokud je matice ryze sloupcově diagonálně dominantní pak Jacobiova iterační metoda konverguje pro každou počáteční aproximaci \(x^0 \in \mathbb{R}^n\)

\begin{align*} |a_{kk}| > \sum_{i = 1, i \neq k}^{n} |a_{ik}| ,\quad k = 1,...,n$ \end{align*}

Gaussova-Seidelova iterační metoda

Z první rovnice vypočteme \(x1\):
\begin{align*} x_1^{k+1} = \frac{1}{a_{11}}(b1 - a_{12}x_2^k - ... - a_{1n}x_n^k) \end{align*}

Z druhé rovnice vypočteme \(x2\), pro \(x1\) použijeme novou iteraci:
\begin{align*} x_2^{k+1} = \frac{1}{a_{22}}(b2 - a_{21}x_1^{k+1} - a_{23}x_3^k - ... - a_{1n}x_n^k) \end{align*}

Z třetí rovnice vypočteme \(x3\), pro \(x1\) a \(x2\) použijeme novou iteraci:
\begin{align*} x_3^{k+1} = \frac{1}{a_{33}}(b3 - a_{31}x_1^{k+1} - a_{32}x_2^{k+1} - a_{34}x_4^k - ... - a_{1n}x_n^k) \end{align*}

Obecně:
\begin{align*} x_i^{k+1} = \frac{b_i}{a_{ii}} - \sum_{j=1}^{i-1} \frac{a_{ij}}{a_{ii}} x_j^{k+1} -\sum_{j=i+1}^{n} \frac{a_{ij}}{a_{ii}} x_j^{k} \quad,\quad i = 1,...,n \end{align*}

Maticový zápis (viz Jacobiova iterační metoda)
\begin{align*} Ax &= b \\ (D+L+U)x &= b \\ (D+L)x &= -Ux + b \\ \end{align*}

\(a_{ii} \neq 0, i = 1,...,n \implies\) matice \(D + L\) je regulární a

\begin{align*} x &= -(D+L)^{-1}Ux +(D+L)^{-1}b \\ \end{align*}

Položme \(T_G = -(D + L)^{-1}U\), Gaussova-Seidelova iterační metoda je tvaru
\begin{align*} x^{k+1} = T_Gx^k + g ,\quad g = (D+L)^{-1}b \end{align*}

Věta o konvergenci

Posloupnost \(\{x_k\}^\infty_{k = 0}\) generovaná metodou \(x^{k+1} = -(D+L)^{-1}Ux +(D+L)^{-1}b\) konverguje pro každou počáteční aproximaci \(x^0 \in \mathbb{R}^n\) právě tehdy, když \(\varrho(T_G) < 1\).

Platí

• Silné řádkové sumační kriterium
• Silné sloupcově sumační kriterium

relaxační metody

Modifikace Gaussovy–Seidelovy metody, \(\omega\) – relaxční parametr

\begin{align*} x_i^{k+1} = (1-\omega)x_i^k + \frac{\omega}{a_{ii}}\left(b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{k+1} -\sum_{j=i+1}^{n} a_{ij} x_j^{k}\right) \end{align*}

Relaxacní metodu lze maticově zapsat jako:
\begin{align*} x^{k+1} &= (D - \omega L)^{-1} [(1-\omega)D + \omega U]x^k + \omega(D - \omega L)^{-1}b \\ T_\omega &= (D - \omega L)^{-1} [(1-\omega)D + \omega U] \end{align*}

Hodnoty parametru \(\omega\):

• Pro \(0 < \omega < 1\) se iterační metody nazývají metodami dolní relaxace. Tyto metody jsou vhodné v případe, že Gaussova-Seidelova metoda nekonverguje.
• Pro \(\omega = 1\) je relaxační metoda totožná s Gaussovou-Seidelovou metodou.
• Pro \(1 < \omega\) se metody nazývají metodami horní relaxace, nebo časteji SOR metodami (Successive Over-Relaxation). Tyto metody lze užít ke zrychlení konvergence Gaussovy-Seidelovy metody.

Věta (Kahan)

Nechť \(a_{ii} \neq 0, i = 1,...,n\) Pak \(\varrho(T_\omega) \ge |\omega - 1|\)
Důsledek: Má smysl uvažovat jen \(\omega \in (0, 2)\)

Věta (Ostrowski-Reich)

Pro pozitivně definitní matici \(A\) platí \(\varrho(T_\omega) < 1\) pro všechna \(\omega \in (0, 2)\)
Třídiagonální matice: \(A = (a_{ij})_{i,j = 1}^n\), \(a_{ij} = 0\), pro \(|i-j| > 1\)

Věta

Nechť \(A\) je třídiagonální pozitivně definitní matice. Pak \(\varrho(T_G) = \varrho^2(T_J) < 1\) a optimální hodnota relaxačního parametru je dána vztahem:
\begin{align*} \omega = \omega_{opt} = \frac{2}{1 + \sqrt{1-\varrho(T_J)}} \end{align*}
Při této volbě je \(\varrho(T_\omega) = |1 - \omega|\)

LU dekompozice

\(A = L \cdot U\)

• \(L\) je dolní trojúhelníková matice
• \(U\) je horní trojúhelníková matice

Obecně \(P \cdot A = L \cdot U\) kde \(P\) je permutační matice, která vyměňuje řádky matice \(A\)
LU dekompozice je v podstatě modifikovaná verze Gaussovy eliminace
\begin{align*} Ax &= b \\ A = LU \implies LUx &= b \\ y = Ux \implies Ly &= b \\ \end{align*}
Řešíme dva systémy rovnic s trojúhelníkovými maticemi.

QR dekompozice

\(A = Q \cdot U\)

• \(Q\) je ortogonální matice, \(Q^{-1} = Q^T\)
• \(U\) nebo nekdy take \(R\) (viz nazev) je horní trojúhelníková matice

QR dekompozice má lepší numerickou stabilitu díky ortogonální transformaci
(co to znamená ? netuším)

\begin{align*} Ax &= b \\ A = QR \implies Ux &= Q^Tb \\ \end{align*}

Choleského dekompozice

Nechť \(A\) je reálná symetrická pozitivně definitní matice: \(A^T = A\)

• \(R\) je horní trojúhelníková matice
  \begin{align*} A = R^T \cdot R \end{align*}

Numerická diferenciace, diferenciační schémata.

Chceme spočítat aproximaci \(f'(x)\)
Pro dané body \(x_0,...,x_n\) a dané funkční hodnoty \(f_0,...,f_n\) kde \(f_k = f(x_k)\)

Polynomy

Prvně spočítáme polynom který lze následně derivovat.
Předpokládámě že zadanými \(n + 1\) body prochází polynom nejvýše \(n\)-tého stupně.

Lagrangeův polynom

\begin{align*} P(x) = \sum_{i=0}^{n}f_iL_i(x) = \sum_{i=0}^{n}f_i\frac{\displaystyle \prod_{j \neq i}(x-x_j)}{\displaystyle \prod_{j \neq i}(x_i-x_j)} \end{align*}

• \(L_i\) – Lagrangeovy základní (base) polynomy
  \begin{align*} L_i(x) = \frac{\pi_i(x)}{\pi_i(x_i)} =\frac{\displaystyle \prod_{j \neq i}(x-x_j)}{\displaystyle \prod_{j \neq i}(x_i-x_j)} \end{align*}

Příklad:

Výpočet základního polynomu \(L_i\) je \(O(n^2)\)
Výpočet celého interpolačního polynomu \(P\) je \(O(n^3)\)

Pro efektivnější výpočet lze využát Hornerovo schéma

\begin{align*} \omega(x) = \displaystyle \prod_{j=0}^{n}(x-x_j) \quad &: O(n^2) \\ \pi_i(x) = \omega(x):(x - x_i) \quad &: O(n)\\ \pi_i(x) \quad &: O(n)\\ P \quad &: O(n^2) \end{align*}

Numerická diferenciace

Příklad 1:

Příklad 2:

Derivace z Taylorova rozvoje

Richardson extrapolation