MNUM Numerické výpočty
Tomáš Raček
Univerzita Mikuláše Koperníka
23. 2. 2013
Něco na zahřátí
1.
d
(︁
T + 𝛿TR
𝛼T Φs
Rd Tst
)︁
dt
=
d
(︁
𝛿TR
𝛼T Φs
Rd Tst
)︁
dt
+
RT
cp
𝜔
Π
+ FT
2.
x3
log(x + 10) = tan(x − 1)2
3.
x cos x = 1
Co mají tyto rovnice společného?
Nelze je řešit analyticky.
ALADIN
ALADIN – předpovědní model počasí používaný v ČR
rovnice teploty:
d
(︁
T + 𝛿TR
𝛼T Φs
Rd Tst
)︁
dt
=
d
(︁
𝛿TR
𝛼T Φs
Rd Tst
)︁
dt
+
RT
cp
𝜔
Π
+ FT
Citát
God does not care about our mathematical difficulties.
He integrates empirically.
Albert Einstein
Přístupy k řešení
Analytický přístup
založen na symbolických úpravách
krajní případy lze jednoduše studovat
nestačí na řešení většiny reálných problémů
Numerický přístup
založen na číselných aproximacích
často poměrně přímočarý
potřeba vizualizace
mezní případy mohou činit obtíže (0/0)
Aplikace numerických metod
chemie (návrh nových struktur, modelování dějů)
biologie (populační, epidemiologické modely)
fyzika (dynamika kapalin)
klimatické modely, předpovědi počasí
astronomie
ekonomie
strojírenství
a mnoho dalších
Nelineární rovnice
uvažujme pouze funkce jedné proměnné
naším cílem je nalézt kořen, tj. řešení rovnice f (x) = 0
jen velmi malou část příkladů lze řešit analyticky
vzorce pro speciální případy (kvadratická, kubická rovnice,. . . )
Naivní přístup k řešení:
zvolme velikost kroku Δx
v rámci intervalu [a, b], ve kterém hledáme kořen,
vyhodnocujeme f (a), f (a + Δx), f (a + 2Δx), . . . , f (b)
pro malá Δx časově velmi náročné
Separace kořene
Věta:
Nechť f je funkce spojitá na [a; b] a f (a) · f (b) < 0.
Pak ∃c ∈ (a; b) : f (c) = 0.
+
a
b
Metoda půlení intervalu
Bisekce (f , l, r, 𝜖)
1. pokud je dosaženo kritérium zastavení, skonči
2. x = (l + r)/2
3. pokud f (l) · f (x) < 0, pak Bisekce (f , l, x, 𝜖),
jinak Bisekce (f , x, r, 𝜖)
Kritéria zastavení
𝜖 – zvolená přesnost
|l − r| < 𝜖
|f ((l + r)/2) | < 𝜖
Vlastnosti
důsledně zachovává separaci kořene → vždy konvergentní
pomalá
Newtonova metoda I – představení
někdy také metoda tečen
vychází z Taylorova rozvoje zanedbáním členů vyšších řádů
f (x) =
∞∑︁
k=0
f (k)(x0)
k!
(x − x0)k
f (x) = f (x0) + f ′
(x0)(x − x0) + . . . +
f (n)(x0)
n!
(x − x0)n
+ Rn
Vzoreček:
xk+1 = xk −
f(xk)
f′(xk)
Newtonova metoda II – vlastnosti
nutná dobrá počáteční aproximace
kvadratická konvergence → jedna iterace zvyšuje počet
platných číslic zhruba dvakrát
užitečná při dolaďování výsledku
X2
ξ
X1 X0
Newtonova metoda III – aplikace
Reciproká odmocnina y = 1/
√
x
často se vyskytuje v aplikacích (např. norma vektoru: v
‖v‖)
obsahuje výpočetně náročné operace – dělení a odmocninu
její aproximaci lze ale spočítat rychle
Newtonova metoda pro y = 1/
√
x
→ hledáme kořen funkce 1/y2 − x = 0.
yk+1 = yk −
1/y2
k − x
(1/y2
k − x)′
= yk −
y−2
k − x
−2y−3
k
yk+1 =
3y−2
k − x
2y−3
k
=
yk(3 − xy2
k )
2
Newtonova metoda IV – Quake III
float Q_rsqrt(float number)
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long *) &y; // evil floating point bit level hacking
i = 0x5f3759df - (i >> 1); // what the f_ck?
y = *(float *) &i;
y = y * (threehalfs - (x2 * y * y)); // 1st iteration
// y = y * (threehalfs - (x2 * y * y));
// 2nd iteration, this can be removed
return y;
}
Obyčejné diferenciální rovnice 1. řádu
Ve zkratce:
obecně tvaru F(x, y, y′) = 0
řešením je funkce y = f (x)
Příklad:
y′
= 2xy
Řešení:
dy
y
= 2x dx y ̸= 0
ln |y| = x2
+ C1 C1 ∈ R
|y| = ex2+C1
= C2 ex2
C2 ∈ R+
y = C3 ex2
C3 ∈ R ∖ {0}
Funkce y = 0 je také řešením → y = Cex2
, C ∈ R.
S-I-R model I – úvod
jednoduchý epidemiologický model
použitelný pro nemoci jako jsou spalničky nebo příušnice
Susceptible – nikdy nepřišli do kontaktu s nemocí
Infectious – šíří nemoc, zůstávají infekční po nějakou dobu
Recovered – již prodělali nemoc a jsou vůči ní imunní
Susceptible Infectious Recovered
S-I-R model II – rovnice
X(t) – množství lidí v kategorii X v čase t
𝛽 – míra kontaktu
𝜈 – míra uzdravení
S-I-R model v rovnicích:
dS
dt
= −𝛽I(t)S(t)
dI
dt
= 𝛽I(t)S(t) − 𝜈I(t)
dR
dt
= 𝜈I(t)
Eulerova dopředná metoda
jednoduchá metoda pro řešení (systémů) dif. rovnic ve tvaru
dy
dt
= f (t, y(t))
postupně počítáme
yn+1 = yn + Δt ·
dy
dt
⃒
⃒
⃒
⃒
tn
y0
y1
y2
y3
y4
S-I-R model III – výsledky
Integrál I – zadání
Spočtěte následující integrál:
∫︁ 1
0
x16
ex−1
dx
Obecněji:
Jn =
∫︁ 1
0
xn
ex−1
dx
Jn =
[︁
xn
ex−1
]︁1
0
−
∫︁ 1
0
nxn−1
ex−1
dx
Jn = 1 − nJn−1
Zřejmě navíc J0 = 1 − 1/e.
Integrál II – implementace
Implementace:
main() {
float Jn = 1.0 - (1.0 / 2.7182818284590452354);
for(int i = 1; i <= 16; i++) {
Jn = 1 - i * Jn;
printf("J_%02d = %12.3f\n", i, Jn);
}
}
Výsledek:
J_01 = 0.368
.
.
.
J_16 = -191438.344
Něco je špatně.
Integrál II – implementace
Implementace:
main() {
float Jn = 1.0 - (1.0 / 2.7182818284590452354);
for(int i = 1; i <= 16; i++) {
Jn = 1 - i * Jn;
printf("J_%02d = %12.3f\n", i, Jn);
}
}
Výsledek:
J_01 = 0.368
.
.
.
J_16 = -191438.344
Něco je špatně.
Reprezentace čísel I – úvod
Celá čísla – datový typ int
32 bitů
rozsah od −2 147 483 648 do 2 147 483 647
resp. 0 až 4 296 967 295 pro neznaménkový typ
odpovídá binárnímu zápisu čísla
v rámci rozsahu naprosto přesné výpočty, mimo něj nesmysly
Reálná čísla
konečná velikost paměti vs. nekonečně mnoho reálných čísel
nutno zvolit odlišný přístup
Reprezentace čísel II – reálná čísla
Reálná čísla – datový typ float
32 bitů
číslo je uloženo ve tvaru (−1)sign · mantissa · 2exponent
sign exponent (8 bits) fraction (23 bits)
02331
0 0 1 1 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 = 0.15625
30 22 (bit index)
nejmenší kladné zobrazitelné číslo: 1, 18 · 10−38
největší zobrazitelné číslo ≈ 3, 4 · 1038
6 platných číslic
speciální hodnoty: ±∞, NaN
Úskalí konečné reprezentace I – porovnávání
Jednoduchá porovnání:
0.5 + 0.5 == 1.0? True.
0.1 + 0.2 == 0.3? False.
některá čísla nelze vyjádřit přesně
(0.1 → 1.0000000149011611e-1)
výsledek je téměř vždy zatížen chybou předchozího výpočtu
operátor == nemá pro reálná čísla moc smysl
Porovnáváme-li dvě reálná čísla, použijme raději test na
dostatečnou blízkost: | x − y | < 𝜖
Úskalí konečné reprezentace II – nepřesnost sčítání
Spočtěte s přesností na tři platné číslice 1, 23 + 0, 001:
výsledek 1, 23
v rámci přesnosti korektní
Ale co v tomto případě?
1, 23 + 0, 001 + . . . + 0, 001
⏟ ⏞
10
= 1, 23
vs.
1, 23 + (0, 001 + . . . + 0, 001) = 1, 24
Sčítání reálných čísel není asociativní.
→
Sčítejme vždy nejdříve čísla se srovnatelným exponentem.
Úskalí konečné reprezentace III – katastrof. zrušení
Uvažme příklad 1, 23 − 1, 22:
binární reprezentace: 1, 0011101 a 1, 0011100
vstupní hodnoty zatíženy chybami 0, 3 %, resp. 0, 1 %
rozdíl binárně 0, 0000001
rozdíl desítkově 0, 0078125
správný výsledek je 0, 01
relativní chyba 22 %!
Snažme se vyhnout odečítání dvou blízkých čísel. Není-li to možné,
počítejme s potenciálně velkou chybou.
Gaussova eliminační metoda (GEM)
řešení systému lineárních rovnic
princip: eliminace prvků pod diagonálou
problém, pokud je vedoucí prvek v aktuálním řádku (pivot)
blízký nule → GEM není stabilní → v praxi se nepoužívá
Stabilita metody
intuitivně: malá chyba na vstupu neovlivní výrazně výsledek
GEM s částečným výběrem vedoucího prvku
prohození dvou řádků výsledek neovlivní
princip: jako vedoucí prvek volím ten, který je v absolutní
hodnotě největší
stabilní
Integrál III – analýza problému
J16 =
∫︁ 1
0
x16
ex−1
dx ̸≈ −191438.344
Kde se stala chyba?
J0 zatížen chybou na vstupu
v každé iteraci násobíme aktuální hodnotou n
rozvoj rekurentní formule pro J16 obsahuje 16!J0
Potřebujeme jiný způsob výpočtu.
Integrál IV – řešení
Jn−1 =
1 − Jn
n
v každém kroku bude chyba redukována faktorem 1/n
nutno nalézt počáteční hodnotu
∫︁ 1
0
xn
ex−1
dx ≤
∫︁ 1
0
xn
dx =
1
n + 1
pro n = 20 platí J20 ≤ 1/21
zvolme J20 = 0
pak J16 ≈ 0, 0557189
správná hodnota 0, 0557193 . . .
Shrnutí hlavních myšlenek
chceme umět popisovat přírodní zákonitosti a děje, na jejich
základě pak dělat predikce, odvozovat vlastnosti ap.
rovnice, které tímto popisem vzniknou, většinou nelze řešit
analyticky
díky masivnímu nárůstu výpočetních kapacit lze ale často
poměrně přímočaře využít numerických metod
je velmi snadné dostat špatný výsledek, pokud si nedám pozor