Chyby zaokrouhlovací#

Zaokrouhlovací chyby vznikají při ukládání desetinných čísel v počítači a při aritmetických operacích s těmito čísli. Kvůli konečné reprezentaci dochází k mnoha případům a numerickým jevům, které můžou vést k nečekanému zhoršení ale i zlepšení výsledku! Některé zajímavé případy si teď ukážeme.

V pythonu se defautně vypisují čísla zkráceně. Vypsat přesnou hodnotu můžeme následovně:

print(0.1, 0.3)
print(format(0.3, ".55f"))

0.1 0.3
0.2999999999999999888977697537484345957636833190917968750

Chyba reprezentace

## DOPLŇTE ##

print(0.1 + 0.2 == 0.3, 0.1+0.2, 0.3)

False 0.30000000000000004 0.3

Vidíme, že součet dvou čísel již zanesl do výsledku malou chybu. Více se k tomuto můžete dozvědět zde.

Odečítání podobných čísel

Odečítání velmi blízkých čísel není vhodné, jelikož při tom vzniká velká relativní chyba. Dvě blízké čísla budou mít zvelké části stejnou mantisu. Po jejich odečtení dostaneme výsledek na malý počet platných číslic:

Pokud $x\approx y$, bude rozdíl $|x\pm y|\approx 0$ a tím pádem relativní chyba bude vysoká. Rozdílem blízkých čísel zesiluje původní chybu.

## DOPLŇTE ##

def f(x):
    print(np.cos(x), 1 - np.cos(x), x**2)
    return (1 - np.cos(x)) / x**2

print(f(1.2e-8))

0.9999999999999999 1.1102230246251565e-16 1.44e-16
0.7709882115452477

Asociativita

V počítači nemusí platit asiciativita sčítání $(a+b)+c=a+(b+c)$. Ilustrujeme to na následující úloze:

## DOPLŇTE ##

base = 0.9
n = 400

S = 0
for i in range(0,n+1):
    S += base**i
print(S)

S = 0
for i in range(n,-1,-1):
    S += base**i
print(S)

9.999999999999993
10.000000000000004

Hromadění chyb 1

Na následujícím příkladu uvidíme, že malá chyba se může ve výpočtu postupně nahromadit.

## DOPLŇTE ##

x:np.float64 = 0.1 + 0.2 - 0.3
print(x)
for i in range(100):
    x = x + x
print(x)

5.551115123125783e-17
70368744177664.0

Hromadění chyb 2

Protipříklad hromadění chyb. Může se zdát, že k velkým chybám dochází až při mnoha výpočtech za sebou. Následující ukázka prokazuje opak.

## DOPLŇTE ##

e = 2.71828
K = 7
for n in np.logspace(1, 7, num=7, base=10):
    arr = np.array([n], dtype=np.float32)
    x = (1 + 1/arr)**arr
    print(x, np.abs(e - x)/e)

[2.593743] [0.04581463]
[2.7048113] [0.00495488]
[2.7170508] [0.00045223]
[2.718597] [0.00011657]
[2.7219622] [0.00135458]
[2.5952268] [0.04526881]
[3.2939677] [0.21178378]

Vzájemné vyrušení chyb

Někdy může dojít k několika chybám během výpočtu, které se v konečném výsledku vyruší. Toho lze využít pro zpřesnění určitých výpočtů.

np.logspace(-5, -16, num=12, base=10)

array([1.e-05, 1.e-06, 1.e-07, 1.e-08, 1.e-09, 1.e-10, 1.e-11, 1.e-12,
       1.e-13, 1.e-14, 1.e-15, 1.e-16])

## DOPLŇTE ##

# 1., 2.
for i in np.logspace(-5, -15, num=12, base=10):
    print((np.exp(i) - 1)/i, (np.exp(i) - 1)/np.log(np.exp(i)))
    
print((np.exp(i) - 1), i, np.log(np.exp(i)))

1.000005000006965 1.0000050000166667
1.000000616442873 1.000000616423623
1.0000000764311845 1.000000075995558
1.0000000062384284 1.000000009369087
0.9999999753484707 1.0000000011550647
0.9999997446031695 1.0000000001424016
0.9999998870652326 1.000000000017556
0.9999996054137743 1.0000000000021643
0.9998187316137019 1.000000000000267
0.9989646455689632 1.0000000000000329
1.012863778612346 1.000000000000004
1.1102230246251565 1.0000000000000004
1.1102230246251565e-15 1e-15 1.110223024625156e-15

Chyby metody#

V případě kdy řešíme úlohu, která v matematické formulaci vyžaduje k vyřešení nekonečný počet kroků, nedokážeme napsat přesný algoritmus. Při převedení na numerickou úlohu se budeme muset dopustit jisté aproximace, aby bylo možné řešení spočítat v konečném počtu kroků.

Typycky výsledná numerická metoda je jistého řádu, určující závislost chyby metody na počtu kroků.

Big-O notace

Tato notace ukazuje řádovou závislost chyby na jistém rozměru úlohy (velikost kroku numerické metody).

Definice:

přávě tehdy když

Tedy big-O notace vyjadřuje, že chyba v určité limitě ($x$ jdoucí k $a$, většinou $a=0$ nebo $a=\mathrm{\infty }$) nepřekročí hodnotu funkce $g$, až na konstantu $M$.

Uvažujme následující úlohu využívající poznatky z matematické analýzy - Taylorův rozvoj funkce:

• $f(x)=\sum _{n=0}^{\mathrm{\infty }}{\displaystyle \frac{f^{(n)}(a)}{n!}}(x-a{)}^n=\sum _{n=0}^{\alpha -1}{\displaystyle \frac{f^{(n)}(a)}{n!}}(x-a{)}^n+O((x-a{)}^{\alpha })$

• nahrazením $h\equiv x-a$

• $f(a+h)=\sum _{n=0}^{\mathrm{\infty }}{\displaystyle \frac{f^{(n)}(a)}{n!}}{h}^n=\sum _{n=0}^{\alpha -1}{\displaystyle \frac{f^{(n)}(a)}{n!}}{h}^n+O(h^{\alpha })$

• Řád metody $\alpha$: chyba je úměrná $h^{\alpha }$

## DOPLŇTE ##

def exp_taylor(x, N):
    fx = 0
    for n in np.arange(0, N):
        fx += x**n / np.math.factorial(n)
    return fx

plt.figure(figsize=(15,5))
plt.subplot(131)
xx = np.linspace(-1, 1, 100)
fe = np.exp(xx)
plt.plot(xx, fe, 'k-', label=r'$e^x$')

plt.subplot(132)
plt.title('absolutní chyba')
plt.subplot(133)
plt.title('relativní chyba')

for i in range(2,4+1):
    ff = exp_taylor(xx, N=i)
    
    plt.subplot(131)
    plt.plot(xx, ff, label=f'stupeň {i}')
    
    abs_err = np.abs(fe - ff)
    rel_err = abs_err / np.abs(fe)
    
    plt.subplot(132)
    plt.plot(xx, abs_err)
    plt.subplot(133)
    plt.plot(xx, rel_err)
    
plt.subplot(131)
plt.legend();

Vylepšená aproximace $e^x$ na $\mathbb{R}$

Taylorův rozvoj funkce aproximuje danou funkci pouze na okolí zvoleného budu. Na aproximaci celé funkce to nestačí.

Co takhle použít Taylorův rozvoj v různých bodech? To se v případě funkce $e^x$ dostaneme do problémů, protože pro výpočet členů rozvoje v bodě $a\ne 0$ potřebujeme znát hodnotu $e^a$, což se právě snažíme spočítat!

My si zde ukážeme malé vylepšení, které již umožní získat hodnotu exponenciální funkce na celé množině $\mathbb{R}$ s rozumnou přesností. V praxi se používá ještě o několik kroků sofistikovanější algoritmus. Více se můžete dozvědět zde.

Vylepšení spočívá v následujícím triku:

kde $x=n\ln 2+x^{\prime }$. Hodnotu $x$ jsme tedy rozdělili celočíselný násobek $\ln 2$, která je nejblíže hodnotě $x$:

Spočítat $2^n$ je v počítači triviální (n-krát bitový posun 1) a $e^{\prime }$ lze spočítat Taylorem, jelikož $x\prime$ je rozumně blízko 0. Tímto trikem jsme výrazně zlepšili přesnost aproximace exponenciální funkce tím, že jsme redukovali maximální hodnotu ($|x^{\prime }|\le \frac{\ln (2)}{2}$) použitou v Taylorově řadě. Jak si ukážeme v následující ukázce, relativní chyba neporoste pro rostoucí $|x|$:

np.log(2)/2

0.34657359027997264

## DOPLŇTE ##

def exp_aproximace_na_r(x, N):
    n = np.round(x / np.log(2))
    x_zbytek = x - n * np.log(2)
    exp_zbytek = exp_taylor(x_zbytek, N) # exp(x')
    #print(n, x_zbytek, exp_zbytek)
    return 2**n * exp_zbytek

plt.figure(figsize=(15,5))
plt.subplot(131)
xx = np.linspace(-15, 15, 100)
fe = np.exp(xx)
plt.plot(xx, fe, 'k-', label=r'$e^x$')

plt.subplot(132)
plt.title('absolutní chyba')
plt.subplot(133)
plt.title('relativní chyba')

for i in range(1,3+1):
    ff = exp_aproximace_na_r(xx, N=i)
    
    plt.subplot(131)
    plt.plot(xx, ff, label=f'stupeň {i}')
    plt.subplot(132)
    abs_err = np.abs(fe - ff)
    plt.plot(xx, abs_err)
    plt.subplot(133)
    rel_err = abs_err / np.abs(fe)
    plt.plot(xx, rel_err)
    
plt.subplot(131)
plt.legend();

Vztah chyb

Je dobré mít na paměti, že i když teoreticky může s menším krokem jít chyba metody k nule, vždy se budou ve výpočtu vyskytovat chyby zaokrouhlovací! Vzájemný vztah těchto dvou chyb ilustruje následující obrázek:

Chyby vstupních dat#

Vstupní data jsou hodnoty, které vkládáme do počítače. Při jejich uložení již dojde ke zaokrouhlovací chybě, která již může způsobit chybný výsledek.

## DOPLŇTE ##

def kvad_koreny(a,b,c):
    D = b**2 - 4*a*c
    print(D)
    if (D > 0):
        return (-b + np.sqrt(D))/(2*a), (-b + np.sqrt(D))/(2*a)
    elif (D < 0):
        return None
    else:
        return -b/a

#b**2 = 4*a*c
print(kvad_koreny(1/3, np.sqrt(4*7*1/3), 7))

1.7763568394002505e-15
(-4.582575631735567, -4.582575631735567)