%% Második gyakorlat: Lebegőpontos számábrázolás, példák stabil/instabil eljárásokra
%
% 1. feladat: Gépeljük be az 1/3 értéket MATLAB-ban. Használjuk a begépelés után a format long és format parancsokat. A help függvény segítségével értelmezzük a format parancs opcióit és kísérletezzünk azokkal!

1/3

format long 
1/3

format
1/3

% Nézzük meg a help segítségével a különböző format opciókat és kísérletezzünk!

help format

% Kísérletezés:

format short
1/3

format shorte
1/3

format shorteng
1/3

format longe
1/3

%
%
% 2. feladat: Lebegőpontos számok esetén az összeadás és  szorzás nem asszociatív művelet, valamint a disztributív tulajdonság sem áll fenn. Győzödjünk meg erről összeadás esetében, ha
% 𝑎 = 0.4, 𝑏 = − 0.5, 𝑐 = 0.1.
% Formálisan logikai operátor segítségével is ellenőrrizük az állítás mivoltát!

format long
a=0.4;
b=-0.5;
c=0.1;

a+(b+c)
(a+b)+c

% Relációk MATLAB-ban:
% egyenlőség: ==,
% nem egyenlő: ~=.
% Logikai értékek:
% igaz: 1,
% hamis: 0.

1+2 == 2+1
1+2 ~= 2+1
1+2 == 1+1

% Alkalmazva a mi példánkra:

a+(b+c)==(a+b)+c

%
% 3. feladat: A help parancs segítségével értelmezzük az eps, realmin és realmax parancsokat. Mit jelentenek ezek az értékek a gépi számok témakörében?

help eps
help realmin
help realmax

eps
realmin
realmax

%
% Gépi számhalmaz:
% M(a,t,k^-, k^+)= \left\{ \pm a^k \left. \sum_{i=1}^t
%  x_i \; a^{-i} \; \right|\; k^- \leq k \leq k^+, \; 0 \leq x_i \leq a-1,\; i=1, \dots t \right\},
% továbbá fel szoktuk tenni, hogy x_1=1 és a nullát pedig mindig hozzáadjuk a halmazhoz.
% \mathcal{M}_{\infty}: a fenti halmaz legnagyobb eleme (a legnagyobb "ábrázolható" szám).
% \varepsilon_0: a legkisebb pozitív elem (a legkisebb "ábrázolható" pozitív szám).
% \varepsilon_1: az 1 és az őt követő szám távolsága (kis egész számok esetén ez a számhalmaz "felbontása").
% Ha egy szám abszolútértéke kisebb, mint \varepsilon_0, akkor nullának tekinthető, ellenkező esetben a hozzá legközelebbi, fenti halmazban szereplő elemmel azonosítjuk, az x számhoz rendelt, halmaz-beli elem jele fl(x).
%
% 4. feladat: (a) Írjuk fel azM(2,2,-2,1)halmaz elemeit! Mi lesz \varepsilon_0, \mathcal{M}_{\infty}, \varepsilon_1?
%                         Megoldás: M(2,2,-2,1)elemei az alábbi halmazzal adhatók meg:
% M(2,2,-2, 1)= \left\{ \pm 2^k \left. \sum_{i=1}^2
%  x_i \; 2^{-i} \; \right|\; -2 \leq k \leq 1, \; 0 \leq x_i \leq 1,\; i=1, 2 \right\} = \\
% = \left\{\pm 2^k \left( \left. \frac12 + \frac{x_k}{4}\right) \right | -2 \leq k \leq 1, x_i \in \{ 0,1\}\right\}.
%                         melynek elemei ezek alapján:
% 0, \; \pm \frac18, \;\pm \frac{3}{16}, \;\pm \frac14, \; \pm \frac38, \; \pm \frac12, \; \pm \frac34, \; \pm 1, \; \pm \frac32.
%                         Most ezeket a számokat tudjuk ábrázolni, így \mathcal{M}_{\infty}=\frac32 (hiszen ez a legnagyobb lehetséges szám), \varepsilon_0=\frac18 (ez a legkisebb pozitív szám) és \varepsilon_1 = \frac32-1 = \frac12 (mivel az 1 rákövetkezője most a \frac32).
%                     (b) Adjuk meg a fenti halmazban az alábbi számokat:
% \frac12, \; \frac{3}{16}, \; \frac{1}{10}, \; \frac{6}{10}, \; \frac{1}{\pi}.
%                             Megoldás: Ezek alapján \mathrm{fl} \left(\frac12\right) = \frac12 (hiszen ő maga is benne van a halmazban), \mathrm{fl} \left(\frac{3}{16}\right) = \frac{3}{16} (ő is benne van a halmazban), \mathrm{fl}\left( \frac{1}{10}\right) = 0 (mivel kisebb, mint \varepsilon_0=\frac18), \mathrm{fl} \left(\frac{6}{10}\right) = \frac{1}{2}, \mathrm{fl} \left(\frac{1}{\pi}\right) = \mathrm{fl} \left(0.3183\dots\right) = \frac{3}{8} (= 0.375).
%
% 5. feladat: Tekintsük az $a_1 + a_2 + \dots + a_n$ alakú összegek kiszámítására az alábbi két algoritmust:
% Algoritmus 1
%
% s=0;
% for i=1:n
%    s=s+a_i;
% end
% Algoritmus 2
%
% s=0;
% for i=n:-1:1
%    s=s+a_i;
% end
% Mit csinál a két algoritmus? Az eddig tanultak fényében a futtatás előtt mit várunk? Írjunk szkriptet, mely kiszámolja az S_n = \sum_{i=1}^n \frac{1}{i(i+1)} értékét mindkét algoritmussal n=999999 esetén!

% Algoritmus 1
format long
n=999999;
s=0;
for i=1:n
    s=s+1/(i*(i+1));
end
s

% Algoritmus 2
n=999999;
s=0;
for i=n:-1:1
    s=s+1/(i*(i+1));
end
s

% Indoklás: az első esetben a kis számokat az összeadás végén adjuk hozzá, amikor az ezekhez tartozó lebegőpontos számok már messze vannak a valódi értékektől, ebből adódik a hiba. A második esetben a kis számoktól indulunk, így ekkor még ezek összeadódnak pontosan az algoritmus elején, ezért lesz ez pontosabb (kvázi az összeg és az éppen hozzáadott tag mindig hasonló nagyságrendű, nincsen nagy eltérés köztük).
%
% 6. feladat: Tekintsük az alábbi integrált:
% I_n = \int_0^1 \frac{x^n}{10+x} dx
% ahol n \in \mathbb{N}.
% Igazoljuk, hogy I_n \geq 0, I_{n+1} \leq I_n, és ezek alapján I_n  \rightarrow 0, ha n \rightarrow 0.
%                 Megoldás:
%                     A pozitivitás nyilvánvaló (pozitív függvény integrálja pozitív).
%                     A második egyenlőtlenség következik abból, hogy x^n \geq x^{n+1} ha x \in (0,1), hiszen pozitív függvények között fentálló egyenlőtlenség esetén integráljaik között is fent fog állni ugyanez a reláció.
%                     A nullához való konvergencia abból következik, hogy a sorozat pozitív, szigorúan monoton csökkenő és nincs pozitív alsó korlátja (ez utóbbi oka, hogy a függvény (x \in (0,1)esetén) pontonkénti nullához tart). Egy másik lehetséges érvelés a rendőr elv használata, melynek során becsüljük az integrált alulról és felülről is:
% \frac{1}{11} \frac{1}{n+1} \leq \int_0^1 \frac{x^n}{11}\leq \int_0^1 \frac{x^n}{10+x} \leq \int_0^1 \frac{x^n}{10} \leq \frac{1}{10} \frac{1}{n+1}
% Igazoljuk, hogy igaz az alábbi rekurzió: I_n = \frac{1}{n} - 10 I_{n-1}, \quad I_0 = ln(1.1).
%                 Megoldás:
%                     A rekurzió átrendezése alapján azt kell igazolnunk, hogy
% \int_0^1 \frac{x^n + 10 x^{n-1}}{10 + x} dx = \frac1n
%                     Egyszerűsítve a bal oldalon (10+x)-el adódik az állítás.
% Írjunk programot, mely adott n-ig kiszámolja a rekurzió lépéseinek értékeit. Mit tapasztalunk nagy n értékek esetén?
% function [I] = rekurzio(n)
%
% % Az n lesz a bemeneti paraméter.
% % Legyen I az a tömb, ami tartalmazza a rekurzió értékeit.
%
% format longe
%
% % Egy n+1 hosszú tömböt veszünk fel, mert a kezdeti, "nulladik"
% % elemet is el fogjuk tárolni. A vektorokat mint egy oszlopú
% % mátrixokat tudjuk definiálni.
%
% I=zeros(n+1,1);
%
% % Fontos, hogy a MATLAB-ban a lista első elemének indexe 1,
% % nem pedig nulla.
%
% I(1) = log(1.1);
%
% % Maga a rekurziós cikulus:
%
% for i=1:n
%     I(i+1) = 1/i - 10*I(i);
% end
%
% end
% Futtassuk ezt a programot (a gyakorlati fájl végén szerkeszthető formában is szerepel)!

rekurzio(32)

% A különös viselkedés oka, hogy a 10-el való szorzások során a gépi hiba minden egyes lépésben tízszeresére nő, így 14 lépés után már nagyobb lesz, mint maga a valós eredmény (a gépi hiba körülbelül {10}^{-16} nagyságrendű, lásd a korábban vizsgált esp parancsot).
% Írjunk programot stabilrekurzió néven, amely az alábbi, "visszafelé" lejátszott rekurziót valósítja meg:
% I_n:=0, \quad I_{n-1}= \frac{1}{10}\left( \frac1n - I_n\right)
%                 Megoldás: A korábbi átrendezéséből adódik ez a szabály. Maga a kód (szerkesztehtő formában a fájl végén):
% function [I] = stabilrekurzio(n)
%
% % Az n lesz a bemeneti paraméter.
% % Legyen I az a tömb, ami tartalmazza a rekurzió értékeit.
%
% format longe
%
% % Egy n+1 hosszú tömböt veszünk fel, mert a kezdeti, "nulladik"
% % elemet is el fogjuk tárolni. A vektorokat mint egy oszlopú
% % mátrixokat tudjuk definiálni.
%
% I=zeros(n+1,1);
%
% % Fontos, hogy a MATLAB-ban a lista első elemének indexe 1,
% % nem pedig nulla.
%
% I(n+1) = 0;
%
% % Maga a rekurziós cikulus:
%
% for i=n:-1:1
%     I(i) = 1/10*(1/i - I(i+1));
% end
%
% end
% Futtassuk a programot!

stabilrekurzio(25)

% Így megkaptuk a stabil rekurziót (hiszen az első elem pont az előző algoritmus első eleme lesz).
%
% 7. feladat: Írjunk programot, mely ábrázolja az 1 környezetében az (x-1)^8 kifejezés értékét az alábbi függvények segítségével!
% f(x)=(x-1)^8
% g(x) =x^8-8x^7+28x^6-56x^5+70x^4-56x^3+28x^2-8x+1
%             Megoldás:

x = linspace(0.98, 1.02, 100);
f=(x-1).^8;
g=x.^8 - 8*x.^7 + 28*x.^ 6 - 56*x.^5 + 70*x.^4 - 56*x.^3 + 28*x.^2 - 8*x + 1;
plot(x, f, 'm*', x, g, 'bo')
legend('f(x)','g(x)')

%             Látható, hogy az összeadások miatt a gépi hibák összeadódnak, ezért nem lesz pontos a g(x) számolás.
%
% 8. feladat: (a) Vizsgáljuk meg az alábbi logikai kifejezések értékeit!
% 2^{66}+1==2^{66},\\
% 2^{66}+10==2^{66},\\
% 2^{66}+100==2^{66},\\
% 2^{66}+1000==2^{66},\\
% 2^{66}+10000==2^{66}.

2^66 + 1 == 2^66
2^66 + 10 == 2^66
2^66 + 100 == 2^66
2^66 + 1000 == 2^66
2^66 + 10000 == 2^66

%                 Látható, hogy csupán a 10000 esetén lesz különbség a kettő értékben, mivel a 2^{66} nagyságrendje 10^{19}, így csupán egy 5 számjegyű szám hozzáadása fogja megváltoztatni (hiszen kb. 15 jegyet tudunk figyelembe venni).
%                     (b) Keressük meg azt a legkisebb n \in \mathbb{N} számot, amelyre a 2^{66} + n == 2^{66} kifejezés értéke hamis!
%                         A logikai értékünk kezdetben igaz, és addig fogjuk növelni az n értékét, amíg hamisat nem kapunk.

n=1;
while 2^66 + n == 2^66
n=n+1;
end
n

%                         Fontos végiggondolni, hogy a while ciklus addig fut, amíg igaz a benne lévő kifejezés: tehát n=8192 esetén még igaz lesz, de n=8193 esetén már nem. Ezt le is ellenőrizhetjük:

2^66 + 8192 == 2^66
2^66 + 8193 == 2^66

%                         Érdemes megjegyezni, hogy nem véletlen, hogy a kapott n=8192 kettőhatvány, hiszen ennél nagyobb számok már befolyásolni fogják a gépi szám utolsó számjegyét.
%
% A 6. feladat programjai:

function [I] = rekurzio(n)

% Az n lesz a bemeneti paraméter.
% Legyen I az a tömb, ami tartalmazza a rekurzió értékeit. 

format longe

% Egy n+1 hosszú tömböt veszünk fel, mert a kezdeti, "nulladik"
% elemet is el fogjuk tárolni. A vektorokat mint egy oszlopú
% mátrixokat tudjuk definiálni.

I=zeros(n+1,1); 

% Fontos, hogy a MATLAB-ban a lista első elemének indexe 1, 
% nem pedig nulla.

I(1) = log(1.1);

% Maga a rekurziós cikulus:

for i=1:n
    I(i+1) = 1/i - 10*I(i);
end

end

function [I] = stabilrekurzio(n)

% Az n lesz a bemeneti paraméter.
% Legyen I az a tömb, ami tartalmazza a rekurzió értékeit. 

format longe

% Egy n+1 hosszú tömböt veszünk fel, mert a kezdeti, "nulladik"
% elemet is el fogjuk tárolni. A vektorokat mint egy oszlopú
% mátrixokat tudjuk definiálni.

I=zeros(n+1,1); 

I(n+1) = 0;

% Maga a rekurziós cikulus:

for i=n:-1:1
    I(i) = 1/10*(1/i - I(i+1));
end

end