Obsah

23. Numerické metody a matematická pravděpodobnost
Diskuze

23. Numerické metody a matematická pravděpodobnost

Numerické rešení algebraických rovnic

Základní pojmy

Mějme soustavu rovnic:

$\matrix{4}{9}{{a_{11}x_1} {+} {a_{12}x_2} {+} {...} {+} {a_{1n}x_n} {=} {b_1} {a_{21}x_1} {+} {a_{22}x_2} {+} {...} {+} {a_{2n}x_n} {=} {b_2} {\vdots} {} {} {} {} {} {} {\vdots}{} {a_{m1}x_1} {+} {a_{m2}x_2} {+} {...} {+} {a_{mn}x_n} {=} {b_m}}$

s neznámými x_1, x_2, ..., x_n

platí Ax = B

metody výpočtu dělíme na přímé a iterační

Příme metody

Cramerovo pravidlo

$x_1 = {D_1}/D, x_2 = {D_2}/D, ..., x_n = {D_n}/D$

… diskriminant

D_i … diskriminant z matice, kde i-tý sloupec nahrazen sloupcem pravých stran

vhodné pro malé množství soustav
existuje právě jedno řešení

příklad

Pomocí Cramerova pravidla najdete rešení soustavy rovnic:

$\matrix{2}{1}{{2x_1 + 3x_2 = 5} {-x_1 + 2x_2 = 8}}$

Řešení:

$D = \delim{|}{\matrix{2}{2}{2 3 {-1} 2}}{|} = 4 + 3 = 7$

$D_1 = \delim{|}{\matrix{2}{2}{5 3 8 2}}{|} = 10 - 24 = -14$

$D_2 = \delim{|}{\matrix{2}{2}{2 5 {-1} 8}}{|} = 16 + 5 = 21$

$x_1 = {D_1}/D = -{14}/7$

$x_2 = {D_2}/D = {21}/7$

Gaussova eliminační metoda

upravujeme matici soustavy na schodovitý tvar a řešení najdeme zpětnout substitucí

Frobeninova věta - soustava rovnic má řešení, pokud hodnot matice soustavy je stejná jako hodnot rozšířené matice soustavy

částečný výběr hlavního prvku - modifikace Gaussovy metody, při úpravě na schodovytý tvar vybíráme jako první nenulový prvek ten, jehož absolutní hodnota je největší (zmenší se tím chyba)

úplný výběr hlavního prvku - podobné částečnému výběru, vybíráme největší hodnoty ve sloupec všech řádků

příklad

Pomocí Gaussovy eliminace vyrešte soustavu rovnic:

$\matrix{3}{1}{ {1.67x_1 − 0.15x_2 + 2.51x_3 = −0.84} {2.15x_1 + 3.02x_2 − 0.17x_3 = 2.32} {1.71x_1 − 2.83x_2 + 1.45x_3 = 1.26} }$

Řešení:

$(\matrix{3}{4}{ {1.67} {−0.15} {2.51} {−0.84} {2.15} {3.02} {−0.17} {2.32} {1.71} {−2.83} {1.45} {1.26} })$

$(\matrix{3}{4}{ {1.67} {−0.15} {2.51} {−0.84} {0} {3.21311} {−3.40144} {3.40144} {0} {−2.67641} {−1.12012} {2.12012} })$

$(\matrix{3}{4}{ {1.67} {−0.15} {2.51} {−0.84} {0} {3.21311} {−3.40144} {3.40144} {0} {0} {-3.95339} {4.95339} })$

$\matrix{3}{7}{ {1.67x1} {−} {0.15x2} {+} {2.51x3} {=} {−0.84} {} {} {3.21311x2} {−} {3.40144x3} {=} {3.40144} {} {} {} {} {−3.95339x3} {=} {4.95339} }$

$x_3 = {4.95339}/{−3.95339} \approx −1.25295$
$x_2 = 1/{3.21311} (3.40144 + 3.40144(−1.25295)) \approx −.26777$
$x_1 = 1/{1.67} (−0.84 + 0.15(−0.26777) − 2.51(−1.25295)) \approx 1.35613$

Iterační metody

Jacobiho metoda

soustavu rovnic převedeme na tvar:

$\matrix{4}{1}{{x_1 = 1/{a_{11}}(b_1 - a_{12}x_2 - a_{13}x_3 - ... - a_{1n}x_n)} {x_2 = 1/{a_{22}}(b_1 - a_{21}x_1 - a_{23}x_3 - ... - a_{2n}x_n)} {\vdots} {{x_n = 1/{a_{nn}}(b_n - a_{n1}x_1 - a_{n2}x_2 - ... - a_{n{n-1}}x_{n-1})}}}$

zvolíme libovolně počáteční řešení:

$x^{(0)} = (x{_1}^{(0)}, x{_2}^{(0)}, ..., x{_n}^{(0)})^T$

dosadíme do rovnic a dostaneme řešení

$x^{(1)} = (x{_1}^{(1)}, x{_2}^{(1)}, ..., x{_n}^{(1)})^T$

obecně:

$\matrix{4}{1}{ {x{_1}^{(r + 1)} = 1/{a_{11}}(b_1 - a_{12}x{_2}^{(r)} - a_{13}x{_3}^{(r)} - ... - a_{1n}x{_n}^{(r)})} {x{_2}^{(r + 1)} = 1/{a_{22}}(b_1 - a_{21}x{_1}^{(r)} - a_{23}x{_3}^{(r)} - ... - a_{2n}x{_n}^{(r)})} {\vdots} {x{_n}^{(r + 1)} = 1/{a_{nn}}(b_n - a_{n1}x{_1}^{(r)} - a_{n2}x{_2}^{(r)} - ... - a_{n{n-1}}x{_{n-1}}^{(r)})} }$

počítá se, dokud řešení nedosáhne např. dané přesnosti nebo není překročen maximální počet kroků

Definice

Matice A se nazývá

řádkově ostře diagonálně dominantní právě tehdy, když

$|a_{ii}| > \sum{j = 1, j \ne i}{n}{|a_{ij}|}$ , pro i = 1 .. n

sloupcově ostře diagonálně dominantní právě tehdy, když

$|a_{jj}| > \sum{i = 1, i \ne j}{n}{|a_{ij}|}$ , pro j = 1 .. n

Věta

Je-li matice A ostře řádkově nebo sloupcove diagonálně dominantní, Jacobiho metoda konverguje

Příklad

Jacobiho metodou rešte soustavu

$\matrix{3}{7}{ {15x_1} {−} {x_2} {+} {2x_3} {=} {30} {2x_1} {−} {10x_2} {+} {x_3} {=} {23} {x_1} {+} {3x_2} {+} {18x_3} {=} {−22} }$

Řešení:

Matice soustavy je diagonálne dominantní, protože platí

$\delim{|}{15}{|} > \delim{|}{−1}{|} + \delim{|}{2}{|}, \delim{|}{−10}{|} > \delim{|}{2}{|} + \delim{|}{1}{|}, \delim{|}{18}{|} > \delim{|}{1}{|} + \delim{|}{3}{|}$

$\matrix{3}{3}{ {{x_1}^{(r+1)}} {=} {1/{15} (30 + {x_2}^{(r)} − 2 {x_3}^{(r)})} {{x_2}^{(r+1)}} {=} {−{1}/{10} (23 − 2{x_1}^{(r)} − {x_3}^{(r)})} {{x_3}^{(r+1)}} {=} {1/{18} (−22 − {x_1}^{(r)} − 3{x_2}^{(r)})} }$

Počáteční aproximace: $\vec{x} = (0, 0, 0)^T$

	${x_1}^{(r)}$	${x_2}^{(r)}$	${x_3}^{(r)}$

Gauss-Seidelova metoda

Podobná Jacobiho metodě, liší se v tom, že při výpočtu další aproximace používá nejnovějších hodnot

$\matrix{4}{1}{ {x{_1}^{(r + 1)} = 1/{a_{11}}(b_1 - a_{12}x{_2}^{(r)} - a_{13}x{_3}^{(r)} - ... - a_{1n}x{_n}^{(r)})} {x{_2}^{(r + 1)} = 1/{a_{22}}(b_1 - a_{21}x{_1}^{(r + 1)} - a_{23}x{_3}^{(r)} - ... - a_{2n}x{_n}^{(r)})} {\vdots} {x{_n}^{(r + 1)} = 1/{a_{nn}}(b_n - a_{n1}x{_1}^{(r + 1)} - a_{n2}x{_2}^{(r + 1)} - ... - a_{n{n-1}}x{_{n-1}}^{(r + 1)})} }$

Definice

Matice se nazývá pozitivně definitní, jestliže pro každý nenulový sloupcový vektor x = (x_1, x_2, ..., x_n)^T platí

x^T A x > 0

Věta

Je-li matice A pozitivně definitní, Gauss-Seidelova metoda konverguje

Příklad

Gauss-Seidelovou metodou rešte soustavu

$\matrix{3}{7}{ {15x_1} {−} {x_2} {+} {2x_3} {=} {30} {2x_1} {−} {10x_2} {+} {x_3} {=} {23} {x_1} {+} {3x_2} {+} {18x_3} {=} {−22} }$

Řešení:

Matice soustavy je diagonálne dominantní, protože platí

$\delim{|}{15}{|} > \delim{|}{−1}{|} + \delim{|}{2}{|}, \delim{|}{−10}{|} > \delim{|}{2}{|} + \delim{|}{1}{|}, \delim{|}{18}{|} > \delim{|}{1}{|} + \delim{|}{3}{|}$

$\matrix{3}{3}{ {{x_1}^{(r+1)}} {=} {1/{15} (30 + {x_2}^{(r)} − 2 {x_3}^{(r)})} {{x_2}^{(r+1)}} {=} {−{1}/{10} (23 − 2{x_1}^{(r+1)} − {x_3}^{(r)})} {{x_3}^{(r+1)}} {=} {1/{18} (−22 − {x_1}^{(r+1)} − 3{x_2}^{(r+1)})} }$

Počáteční aproximace: $\vec{x} = (0, 0, 0)^T$

	${x_1}^{(r)}$	${x_2}^{(r)}$	${x_3}^{(r)}$

Numerické metody rešení nelineárních rovnic

Řešením rozumíme každé

vyhovující rovnici f(x) = 0

Věta

Je-li funkce spojitá na intervalu <a, b> a platí

f(a).f(b) < 0 ,

pak v intervalu <a, b> leží alespoň jeden kořen rovnice f(x) = 0 .

Metoda půlení intervalů

Mějme fci

a interval <a, b>

, na kterém platí f(a).f(b) < 0

1. Zvolme střed intervalu $x_0 = {a + b}/2$ . Dostaneme intervaly ,
2. Vyberme ten interval, na kterém je opět splněna podmínka , označme jej
3. vraťmě se ke kroku 2

Skončíme, pokud platí: $b_k - a_k < 2\epsilon$

Řešením je potom $x_k = {b_k + a_k}/2$

Příklad

Metodou pulení intervalu najdete kladný koren rovnice

e^x + x^2 − 3 = 0

s přesností $\epsilon = 0.01$ .

Řešení:

Zvolme interval <0; 1>

b_6 − a_6 < 2 0.01

$x_6 \approx 0.84$

Metoda Regula Falsi

Myšlenka - sestrojíme úsečku s krajními body f(a), f(b)

, určíme průsečík s osou x. Vybereme interval, ve kterém leží kořen rovnice. Opět sestrojíme úsečku v krajních bodech. Opět určíme průsečík a interval, ve kterém je řešení.

Podobně jako metoda půlení intervalů.

Mějme fci a interval <a, b> , na kterém platí f(a).f(b) < 0 .

1. spočtěme přiblížné řešení podle vzorce: $x_k = b_k - {b_k - a_k}/{f(b_k) - f(a_k)} f(b_k)$
2. dostame dva intervaly , , vyberme ten, pro nějž je splněna podmínka
3. vraťmě se na krok 1.

Počítáme, dokud neplatí: $x_{k+1} - x_{k} < \epsilon$

Příklad

Metodou regula falsi najdete kladný koren rovnice

e^x + x^2 − 3 = 0

s přesností $\epsilon = 0.01$ .

Řešení:

Zvolme interval <0; 1>

$\delim{|}{x_2 − x_1}{|} < 0.01$

$x_2 \approx 0.83$

Newtonova metoda (metoda tečen)

Myšlenka - sestrojíme tečnu v krajním bodu intervalu (ten s větší fčí hodnotou) a určíme průsečík s osou x. V tomto bodě vypočteme fčí hodnotu a z něj sestrojíme tečnu. Celý proces se opakuje.

Mějme fci a interval <a, b> , na kterém platí f(a).f(b) < 0 .

průsečík s osou spočítáme podle vztahu: $x_{k + 1} = x_k - {f(x_k)}/{f\prime (x_k)}$

Počítáme, dokud neplatí: $x_{k+1} - x_{k} < \epsilon$

Věta (Fourierova podmínka)

Nechť v intervalu <a, b> leží jediný kořen rovnice f(x) = 0 a nechť $f\prime (x)$ a $f\prime\prime (x)$ jsou spojité a nemění znaménko na intervalu <a, b> . Zvolíme-li za počáteční aproximaci $x_0 \in <a, b>$ tak, aby byla splněna podmínka

$f(x_0)f\prime\prime (x_0) > 0$ ,

Newtonova metoda bude konvergovat.

Příklad

Newtonovou metodou najdete záporný koren rovnice

e^x + x^2 − 3 = 0

s přesností $\epsilon = 0.01$ .

Řešení:

Zvolme interval <-2; -1>

f(x) = e^x + x^2 − 3

$f\prime(x) = e^x + 2x$

$f\prime\prime(x) = e^x + 2$

Na celém intervalu <−2, −1> je $f\prime(x) < 0$ a $f\prime\prime(x) > 0$

Protože $f(−2) = e^{−2} + 1 > 0$ a $f(−1) = e^{−1} − 2 < 0$ , zvolíme x_0 = −2

$x_{k + 1} = x_k − {f(x_k)}/{f\prime(x_k)} = x_k − {e^{x_k} + {x_k}^2 − 3}/{e^{x_k} + 2x_k}$

$\matrix{4}{1}{ {x_0 = −2} {x_1 \approx −1.70623} {x_2 \approx −1.67752} {x_3 \approx −1.67723} }$

$\delim{|}{x_3 − x_2}{|} < 0.01$

$x_3 \approx −1.68$

Metoda prosté iterace

vhodné pro rovnice, které lze vyjádřit ve tvaru x = g(x)

Zvolíme počáteční bod aproximace x_0 , x_k potom spočteme podle vztahu $x_k = g(x_{k - 1})$

Počítáme, dokud neplatí: $x_{k+1} - x_{k} < \epsilon$

Věta

Nechť funkce g zobrazuje interval <a, b> do sebe a má na tomto intervalu derivaci. Jestliže existuje číslo $\alpha \in <0, 1)$ tak, že

$|g\prime (x)| < \alpha, \forall x \in <a, b>$ ,

potom iterační metoda konverguje

Příklad

Metodou prosté iterace najdete záporný koren rovnice

e^x + x^2 − 3 = 0

s přesností $\epsilon = 0.01$ .

Řešení:

Zvolme interval <-2; -1>

$x_2 = 3 − e^x \doubleright x = ±\sqrt{3 − e^x}$

$g(x) = −\sqrt{3 − e^x}$

Ověření podmínky konvergence:

$g\prime(x) = {e^x}/{2\sqrt{3 − e^x}}$

maximem $g\prime(x)$ je $g\prime(-1) \le 0.12 < 1$ ,

a protože

$g(−2) \approx −1.69 \in <−2,−1>$ a $g(−1) \approx −1.62 \in <−2,−1>$ ,

metoda konverguje.

Zvolme x_0 = −2

$x_{k + 1} = g(x_k) = −\sqrt{3 − e^{x_k}}$

$\matrix{4}{3}{ {x_0} {=} {−2} {x_1} {\approx} {−1.69253} {x_2} {\approx} {−1.67808} {x_3} {\approx} {−1.67728} }$

$\delim{|}{x_3 − x_2}{|} < 0.01$

$x_3 \approx −1.68$

Numerické metody rešení soustav nelineárních rovnic

Nechť máme:

$\matrix{4}{3}{{f_1(x_1, x_2, ..., x_n)} {=} {0} {f_2(x_1, x_2, ..., x_n)} {=} {0} {} {\vdots} {} {f_n(x_1, x_2, ..., x_n)} {=} {0}}$

$\vec{F(\vec{x})} = \vec{0}$

$\vec{F} = (f_1, ..., f_n)^T$ , $\vec{x} = (x_1, ..., x_n)^T$

Metoda prosté iterace

Soustavu rovnic upravíme na tvar:

$\matrix{4}{1}{ {x_1 = g_1(x_1, x_2, ..., x_n)} {x_2 = g_2(x_1, x_2, ..., x_n)} {\vdots} {x_n = g_n(x_1, x_2, ..., x_n)} }$

$\vec{x} = \vec{G(vec{x})}$

$\vec{G} = (g_1, ..., g_n)^T$

Analogie u jedné rovnice, zvolíme počáteční řešení $\vec{x^(0)}$ a počítáme:

$\vec{x^{(k + 1)}} = \vec{G(\vec{x^{(k)}})}$

Počítá se dokud neplatí: $|| \vec{x^{(k)}} − \vec{x^{(k−1)}}|| < \epsilon$

Najít vhodné iterační funkce muže být velmi obtížné. Proto se daleko častěji používá Newtonova metoda.

Dále označme

$\vec{G\prime} = (\matrix{4}{4}{ {{\partial g_1}/{\partial x_1}} {{\partial g_1}/{\partial x_2}} {...} {{\partial g_1}/{\partial x_n}} {{\partial g_2}/{\partial x_1}} {{\partial g_2}/{\partial x_2}} {...} {{\partial g_2}/{\partial x_n}} {} {} {\ddots} {} {{\partial g_n}/{\partial x_1}} {{\partial g_n}/{\partial x_2}} {...} {{\partial g_n}/{\partial x_n}} })$

Věta

Nechť zobrazuje uzavřenou oblast D do sebe a je v této oblasti diferencovatelná. Jestliže existuje číslo $\alpha \in <0, 1)$ tak, že

$||G\prime|| < \alpha, \forall x \in D$ ,

kde $||G\prime||$ je řádková nebo sloupcová norma matice $G\prime$ , pak metoda diverguje.

Příklad

Metodou prosté iterace najdete koren soustavy rovnic

$\matrix{2}{3}{ {3_x + x^2y − 3} {=} {0} {x^2 − 5y} {=} {0} }$ ,

který leží v oblasti D = <1/2, 1> * <0, 1/2> s přesností 0,01 .

Řešení:

Zvolme např.

$g_1(x, y) = 1 − {x^2y}/3$ $g_2(x, y) = {x^2}/5$

G = (g_1, g_2)^T zobrazuje D do sebe:

Jestliže $x \in <1/2, 1>$ a $y \in <0, 1/2>$ , pak ${x^2y}/3 \in <0, 1/6>$ a tedy $g_1(x, y) \in <5/6, 1> \subset <1/2, 1>$ .

Podobně $g_2(x, y) \in <1/{20}, 1/5> \subset <0, 1/2>$ .

Konvergence:

Oveříme, zda ${||G\prime||} < 1$ , neboli zda $\delim{|}{{\partial g_1}/{\partial x}}{|} + \delim{|}{{\partial g_1}/{\partial y}}{|} \le \alpha$ i $\delim{|}{{\partial g_2}/{\partial x}}{|} + \delim{|}{{\partial g_2}/{\partial y}}{|} \le \alpha$ .

$G\prime = (\matrix{2}{2}{ {-{2xy}/3} {−{x^2}/3} {{2x}/5} {0} })$

Jestliže $x \in <1/2, 1>$ a $y \in <0, 1/2>$ , pak $\delim{|}{−{2xy}/3}{|} + \delim{|}{−{x^2}/3}{|} \le 1/3 + 1/3 = 2/3 < 1$ a $\delim{|}{{2x}/5}{|} \le 2/5 < 1$ . (Tedy $\alpha = 2/3$ .) Podmínky konvergence jsou splneny.

Jako pocátecní aproximaci mužeme zvolit napr. (x_0, y_0) = (1, 0) .

$x_{k + 1} = g_1(x_k, y_k) = 1 − {{x_k}^2{y_k}}/3$
$y_{k + 1} = g_2(x_k, y_k) = {{x_k}^2}/5$


0	1	0
1	1	0.2
2	0.933	0.2
3	0.942	0.174
4	0.948	0.177

$\delim{|}{x_4 − x_3}{|} < 0.01$ i $\delim{|}{y_4 − y_3}{|} < 0.01 \doubleright (x, y) \approx (0.95, 0.18)$

Newtona metoda

$\vec{F} = (f_1, ..., f_n)^T$

Vychází se ze vztahu:

$\vec{x^{(k+1)}} = \vec{x^{(k)}} − (\vec{F\prime (\vec{x^{(k)}}))}^{-1} \vec{F(\vec{x^{(k)}})}$

Je potřeba spočítat inverzní matici, což je pracné, proto vztah upravíme:

$\vec{F\prime (\vec{x^{(k)}})}(\vec{x^{(k+1)}} - \vec{x^{(k)}}) = -\vec{F(\vec{x^{(k)}})}$

Označme ${\delta}^{(k)} = x^{(k+1)} − x^{(k)} = ({\delta_1}^{(k)}, ..., {\delta_n}^{(k)})^T$

Řešíme potom soustavu:

$\vec{F\prime (\vec{x^{(k)}})}{\delta}^{(k)} = -\vec{F(\vec{x^{(k)}})}$

s neznámými ${\delta_1}^{(k)}, ..., {\delta_n}^{(k)}$ .

Novou aproximaci potom dostaneme:

$x^{(k+1)} = x^{(k)} + {\delta}^{(k)}$

Počítáme, dokud není splněno:

$||x^{(k)} − x^{(k−1)}|| < \epsilon$ neboli $||{\delta}^{(k−1)}|| < \epsilon||$

Příklad

Newtonovou metodou najdete rešení soustavy rovnic

$\matrix{2}{3}{ {(x − 1)^2 + y^2 − 4} {=} {0} {x + (y + 1)^2 − 1} {=} {0} }$

s presností $\epsilon = 0.01$ .

Řešení:

Zvolme x^(0) = (0,−2) .

f_1(x, y) = (x − 1)^2 + y^2 − 4
f_2(x, y) = x + (y + 1)^2 − 1

$F\prime = (\matrix{2}{2}{ {{\partial f_1}/{\partial x}} {{\partial f_1}/{\partial y}} {{\partial f_2}/{\partial x}} {{\partial f_2}/{\partial y}} }) = (\matrix{2}{2}{ {2(x − 1)} {2y} {1} {2(y + 1)} })$

1. krok

$F\prime(0,−2) = (\matrix{2}{2}{ {−2} {−4} {1} {−2} })$

$F(0,−2) = (\matrix{2}{1}{1 0})$

$\matrix{2}{5}{ {−2\delta_1} {−} {4\delta_2} {=} {−1} {\delta_1} {−} {2\delta_2} {=} {0} }$

$\delta_1 = 1/4 = 0.25, \delta_2 = 1/8 = 0.125$

x^(1) = (0 + 0.25,−2 + 0.125) = (0.25, −1.875) .

2. krok

$\delta_1 = 0.01225$ , $\delta_2 = 0.01593$

x^(2) = (0.26225, −1.85907) .

3. krok

$\delta_1 = −0.00004$ , $\delta_2 = 0.00012$

x^(3) = (0.26221, −1.85894) .

$\delim{|}{\delta_1}{|} < 0.01$ a $\delim{|}{\delta_2}{|} < 0.01$ , můžeme skončit.

Numerické rešení obyčejných diferenciálních rovnic

Numerické derivování a integrování

derivování

Má-li funkce

druhou derivaci na intervalu <x_0, x_1>

, pak existují body $\xi_0, \xi_1 \in <x_0, x_1>$ tak, že platí

$f\prime(x_0) = {f(x_1) - f(x_0)}/h - h/2 f\prime\prime (\xi_0)$

$f\prime(x_1) = {f(x_1) - f(x_0)}/h - h/2 f\prime\prime (\xi_0)$

Numerické derivování - vzorce

integrování

Integrály funkce jedné proměnné

Numerické rešení diferenciálních rovnic

Pojmy, předpoklady

Máme interval <a, b>

, který rozdělíme na

dílků: $a = x_0 < x_1 < ... < x_{n - 1} = b$ , x_i

nazýváme uzlové body

$h_i = x_{i + 1} - x_i$ … krok

Podle toho, jestli metoda využívá předchozí hodnoty nebo ne, dělíme metody na:

jednokrokové
vícekrokové ( $y_{n-1}, y_{n-2}, ...$ )

Počáteční úloha

$y\prime(x) = f(x, y), y(x_0) = y_0$

Hledáme fci y(x)

Eulerova metoda

Rovnici $y\prime(x) = f(x, y), y(x_0) = y_0$ upravíme:

${y(x_{i + 1}) - y(x_i)}/h \approx f(x, y), y(x_0) = y_0$ .

Další úpravou získáme:

$y_{i + 1} = y_i + hf(x_i, y_i)$

Příklad

Eulerovou metodou s krokem h = 0.1 řešte počáteční úlohu

$y\prime = x^2 − y, y(0) = 1$

na intervalu <0, 0.5> .

Řešení:

$y_{i+1} = y_i + 0.1({x_i}^2 − y_i)$

0	1	2	3	4	5
0	0.1	0.2	0.3	0.4	0.5
1	0.9	0.811	0.7339	0.6695	0.6186
1	0.9052	0.8213	0.7492	0.6897	0.6435

Modifikace eulerovy metody

První modifikace:

k_1 = f(x_n, y_n)

k_2 = f(x_n + 1/2 h, y_n + 1/2 h k_1)

$y_{n + 1} = y_n + h k_2$

Druhá modifikace:

k_1 = f(x_n, y_n)
k_2 = f(x_n + h , y_n + h k_1)
$y_{n + 1} = y_n + 1/2 h(k_1 + k_2)$

Runge-Kuttovy metody

Obecný tvar:

$y_{n + 1} = y_n + h(w_1 k_1 + ... + w_s k_s)$ ,

kde

k_1 = f(x_n, y_n)

$k_i = f(x_n + \alpha_i h, y_n + h \sum{j = 1}{i - 1}{\beta_{ij}k_j}), i = 2, ..., s$

Nejznámější Runge-Kuttova metoda 4. řádu:

$y_{n + 1} = y_n + 1/6 h(k_1 + 2k_2 + 2k_3 + k_4)$
k_1 = f(x_n, y_n)
k_2 = f(x_n + 1/2 h, y_n + 1/2 hk_1)
k_3 = f(x_n + 1/2 h, y_n + 1/2 hk_2)
k_4 = f(x_n + h, y_n + hk_3)

Příklad

Rungovou-Kuttovou metodou řešte počáteční úlohu

$y\prime = x^2 − y, y(0) = 1$

s krokem h = 0.1 na intervalu <0, 0.5> .

Řešení:

$\matrix{5}{3}{ {k_1} {=} {f(0, 1) = 0^2 − 1 = −1} {k_2} {=} {f(0 + 1/2 0.1, 1 + 1/2 0.1(−1)) = f(0.05, 0.95) = −0.9475} {k_3} {=} {f(0 + 1/2 0.1, 1 + 1/2 0.1(−0.9475)) = f(0.05, 0.952625) = −0.950125} {k_4} {=} {f(0 + 0.1, 1 + 0.1(−0.950125)) = f(0.1, 0.9049875) = −0.8949875} {y_1} {=} {y_0 + 1/6 0.1(k_1 + 2k_2 + 2k_3 + k_4) \approx 0.9051627} }$


0	0	1	1	$\matrix{4}{1}{0 0.05 0.05 0.1}$	$\matrix{4}{1}{1 0.95 0.952625 0.9049875}$	$\matrix{4}{1}{{k_1 = −1} {k_2 = −0.9475} {k_3 = −0.950125} {k_4 = −0.8949875}}$
1	0.1	0.9051627	0.9051626	$\matrix{4}{1}{0.1 0.15 0.15 0.2}$	$\matrix{4}{1}{0.9051627 0.8604046 0.8632675 0.8210860}$	$\matrix{4}{1}{{k_1 = −0.8951627} {k_2 = −0.8379046} {k_3 = −0.8407675} {k_4 = −0.7810860}}$
2	0.2	0.8212695	0.8212693	$\matrix{4}{1}{0.2 0.25 0.25 0.3}$	$\matrix{4}{1}{0.8212695 0.7822060 0.7852842 0.7489911}$	$\matrix{4}{1}{{k_1 = −0.7812695} {k_2 = −0.7197060} {k_3 = −0.7227842} {k_4 = −0.6589911}}$
3	0.3	0.7491822	0.7491818	$\matrix{4}{1}{0.3 0.35 0.35 0.4}$	$\matrix{4}{1}{0.7491822 0.7162230 0.7194960 0.6894826}$	$\matrix{4}{1}{{k_1 = −0.6591822} {k_2 = −0.5937230} {k_3 = −0.5969960} {k_4 = −0.5294826}}$
4	0.4	0.6896804	0.6896800	$\matrix{4}{1}{0.4 0.45 0.45 0.5}$	$\matrix{4}{1}{0.6896804 0.6631964 0.6666456 0.6432659}$	$\matrix{4}{1}{{k_1 = −0.5296804} {k_2 = −0.4606964} {k_3 = −0.4641456} {k_4 = −0.3932659}}$
5	0.5	0.6434699	0.6434693

Rešení soustav diferenciálních rovnic

$\matrix{4}{2}{ {y_1\prime = f_1(x, y_1, y_2, ..., y_n)} {y_1(x_0) = \eta_1} {y_2\prime = f_2(x, y_1, y_2, ..., y_n)} {y_2(x_0) = \eta_2} {\vdots} {\vdots} {y_n\prime = f_1(x, y_1, y_2, ..., y_n)} {y_n(x_0) = \eta_n} }$

$\vec{y\prime} = \vec{f(\vec{x}, \vec{y})}, \vec{y(\vec{x_0})} = \vec{\eta}$

$\vec{y} = (y_1, ..., y_n)^T, \vec{f} = (f_1, ..., f_n)^T, \vec{\eta} = (\eta_1, ..., \eta_n)^T$

Eulerova metoda

$\vec{y_{n+1}} = \vec{y_n} + h\vec{f(x_n, \vec{y_n})}$

Runge-Kutte 4. řádu

$\vec{y_{n + 1}} = \vec{y_n} + 1/6 h(\vec{k_1} + 2\vec{k_2} + 2\vec{k_3} + \vec{k_4})$
$\vec{k_1} = \vec{f(x_n, \vec{y_n})}$
$\vec{k_2} = \vec{f(x_n + 1/2 h, \vec{y_n} + 1/2 h\vec{k_1})}$
$\vec{k_3} = \vec{f(x_n + 1/2 h, \vec{y_n} + 1/2 h\vec{k_2})}$
$\vec{k_4} = \vec{f(x_n + h, \vec{y_n} + h\vec{k_3})}$

Rešení diferenciálních rovnic vyššího rádu

ODE

$y^{(n)} = f(x, y, y\prime, . . . , y^{(n−1)}), y(x_0) = y_0, y\prime(x_0) = {y\prime}_0, ..., y^{(n−1)}(x_0) = {y_0}^{(n−1)}$

mužeme převést na soustavu diferenciálních rovnic prvního řádu:

Označme:

$\matrix{4}{2}{ {{y\prime}_1 = y_2} {y_1(x_0) = y_0} {{y\prime}_2 = y_3} {y_2(x_0) = {y\prime}_0} {\vdots} {\vdots} {{y\prime}_n = f(x, y_1, y_2, ..., y_n)} {y_n(x_0) = {y^{(n - 1)}}_0} }$

Rozložení pravděpodobnosti

Prerekvizity

Pravděpodobnost
Diskrétní náhodná veličina
Spojitá náhodná veličina
Střední hodnota
Rozptyl

Diskrétní rozložení pravděpodobnosti

Binomické rozložení

Pravděpodobnost úspěchu je

, pravděpodobnost neúspěchu 1 − p

. Náhodná veličina

, která udává počet výskytů úspěchu při

nezávislých opakování experimentů, má tzv. binomické rozdělení pravděpodobnosti.

$P(X = r) = (\matrix{2}{1}{{N} {r}}) p^r (1 - p)^{N - r}$

Zapsano: $X \approx Bi(N, p)$ .

E(X) = Np

D(X) = Np(1 - p)

Příklad

Stanovte pravděpodobnost jevu, že ve 4 hodech kostkou padnou právě 3 šestky.

Řešení:

$P(X = i) = (\matrix{2}{1}{n i}) p^i (1 - p)^{n-i}$ ,

$P(X = 3) = (\matrix{2}{1}{4 3}) (1/6)^i (1 - 1/6)^{n-i} \approx 0.015$

Poissonovo rozdelení

Rozložení nazýváme poissonovské, pokid platí:

Pravděpodobnost výskytu události v intervalu závisí pouze na , nikoli na počtu událostí, které nastaly před okamžikem , ani na samotném
Pravděpodobnost, že v časovém intervalu délky k výskytu žádné události nedojde, je kladná, ale menší než 1.
Pro malá nastane v intervalu délky nejvýše jedna událost

$P(Y = k) = {{\lambda}^k}/{k!} e^{-\lambda}$ pro k = 0, 1, 2, ...

$E(Y) = D(Y) = \lambda$

Příklad

Na určítém území dopadí každoročně v určitém období průměrně 7 meteoritů. Jaká je pravděpodobnost, že daný rok v tomto období dopadnou 3 meteority.

Řešení:

$P(Y = 3) = {7^3 e^{-7}}/{3!} \approx 0.052129$

Spojité náhodné veličiny

Rovnomerné rozdelení

$f(t) = \lbrace \matrix{2}{1}{{1/{b - a} ... t \in <a, b>} {0 ... jinak}}$

$F(t) = \lbrace \matrix{3}{1}{{0 ... t \le a} {{t - a}/{b - a} ... t \in (a, b)} {1 ... t \ge b}}$

E(X) = {a + b}/2

$D(X) = {(b - a)^2}/{12}$

Ro(od, do)

Normální rozdelení

$f(t) = 1/{\sqrt{2\pi}\sigma} e^{-{(t - \mu)^2}/{2{\sigma}^2}}$

$E(X) = \mu$

$D(X) = {\sigma}^2$

$N(\mu, \sigma^2)$

U-rozložení

Pro počítání pravděpodobnosti normální rozdělení pro P(x_1 < X < x_2)

Distribuční funkci nelze zintegrovat, proto používáme převod na U-rozložení.

$U = {X - \mu_x}/{\sigma_x}$

Potom platí:

$\phi(u) = P(U < u) = \int{-\infty}{u}{1/{\sqrt{2\pi}} e^{-{t^2}/{2}} dt}$

$\phi(-u) = 1 - \phi(u)$

Exponenciální rozložení

$f(x) = \lambda e^{-\lambda x}, \lambda > 0$ pro x > 0

a nulové pro $x \le 0$

$E(x) = 1/{\lambda}$

$D(x) = 1/{\lambda^2}$

Generování pseudonáhodných čísel

Základem pro generování jakéhokoli rozložení je generátor rovnoměrně rozložených čísel v intervalu <0, 1). Nejčastějí se používají generátory využívající princip lineárního kongruentního generátoru:

$x_{i + 1} = (ax_i + b) mod m$ ,

kde a, b, m jsou vhodně zvolené konstanty.

Generátor generuje čísla v rozsahu $0 \le x_i \le m$ .Pro převod na interval <0, 1) je potřeba výsledné pseudonáhodné číslo podělit , též perioda generátoru.

Pro dobré statistické výsledky je potřeba zvolit vhodné konstanty. Na počítačích se volí m = 2^n , kde je počet bitů typu unsigned integer (operace modulu se potom děje automaticky).

Nevýhody:

zhoršení statistických vlastností generátoru špatnou volbou konstant
závislost po sobě jdoucích generovaných čísel u lineárních generátorů
pro převod na je potřeba vydělit
málo náhodné poslední významné bity

příklad generátoru

static unsigned long ix = seed; // počáteční_hodnota
double Random(void) {
   ix = ix * 69069L + 1; // implicitní operace modulo
   return ix / ((double)ULONG_MAX + 1);
}

Další možností je Mersenne twister. Nemá většinou problémů předchozího generátoru.

Transformace rozložení

Převod rovnoměrného rozložení na jiné, které požadujeme.

metody

inverzní transformace
vylučovací metoda
kompoziční metoda

inverzní transformace

využívá inverzní fuknci k distribuční fci požadovaného rozložení
vhodná tam, kde je snadné inverzní fci spočítat

postup:

1. inverze distribuční funkce: $F^{-1}$
Generování náhodného čísla s rovnoměrným rozložením:
Výsledek: $y = F^{−1}(x)$

Příklad

$F(x) = 1 − e^{−{{x−A}/b}}$

$F^{−1}(x) = A − b ∗ ln (1 − x)$

pro A = 0, b = 1

Generování:

y_i = A − b ∗ ln (1 − x_i)

vylučovací metoda

Máme graf fce hustoty pravděpodobnosti. Náhodně generujeme body s rovnoměrným rozložením v obdélníku daného grafu, pokud bod leží pod grafem, vrátíme ho, jinak opakujeme generování.