%% Ötödik gyakorlat
%% Numerikus deriválás és integrálás
%
%% Numerikus deriválás
% Alapgondolat: Keressük egy f(x) függvény deriváltjának értékét egy adott x = x_0 pontban úgy, hogy f(x) alakja nem ismert, csak az értékei bizonyos pontokban.
% Híres közelítések:
% Jobboldali/haladó differenciahányados:
% f'(x_0) \approx \frac{f(x_0+h)-f(x_0)}{h}
% Baloldali/retrográd differenciahányados:
% f'(x_0) \approx \frac{f(x_0) - f(x_0-h)}{h}
% Centrális/központi differenciahányados:
% f'(x_0) \approx \frac{f(x_0+h) - f(x_0-h)}{2h}
% Centrális közelítés második deriváltra:
% f''(x_0) \approx \frac{f(x_0-h) - 2 f(x_0)  + f(x_0+h)}{h^2}
%
% 1. feladat: Bizonyítsuk be, hogyha f \in C^3 (azaz f háromszor folytonosan deriválható), akkor
% \frac{-3 f(x_0) + 4 f(x_0+h) - f(x_0+2h)}{2h} = f'(x_0) + O(h^2).
% Megoldás: Írjuk fel az x_0 középpontú Taylor sort az x_0+h és x_0+2h pontokban!
% f(x_0+h) = f(x_0) + f'(x_0) h + \frac{f''(x_0)h^2}{2} + O(h^3)
%  f(x_0+2 h) = f(x_0) + f'(x_0) 2 h + \frac{f''(x_0)(2 h)^2}{2} + O(h^3)
% Ezeket behelyettesítve a fenti formulába:
% \frac{1}{2h} \left( - 3 f(x_0) + 4 \left( f(x_0) + f'(x_0) h + \frac{f''(x_0)h^2}{2} + O(h^3)\right) - \right.\\
% \qquad \qquad \qquad\left.- \left( f(x_0) + f'(x_0) 2 h + \frac{f''(x_0)(2 h)^2}{2} + O(h^3)\right)\right) =
% = \frac{1}{2h} \left( f(x_0) (-3+4-1) + f'(x_0) (4h - 2h) + f''(x_0) \left(4 \frac{h^2}{2} - \frac{4 h^2}{2}\right) + O(h^3) \right)=
% = \frac{1}{2h} (2h f'(x_0) + O(h^3)) = f'(x_0) + O(h^2).
%
% 2. feladat: (a) Írjunk MATLAB programot, mely a \sin''(0.5) = -0.479425538604203\dots értéket a másodrendű centrális differenciával közelíti adott h lépésköz mellett! Számítsuk ki a hibát is!
% Megoldás:

h = 0.1;

kozelites = (sin(0.5-h) - 2*sin(0.5) + sin(0.5+h))/h^2
hiba = abs(kozelites-(-0.479425538604203))

%
% (b) Vizsgáljuk meg, mennyi lesz a hiba 10^{-2}, 10^{-3}, 10^{-4}, 10^{-5} és 10^{-6} lépésközök esetén! Mit tapasztalunk? Mivel lehet magyarázni a jelenséget?
% Megoldás:

format long

lepesek = [10^(-2),10^(-3),10^(-4),10^(-5),10^(-6)];
hibak = zeros(5,1);

for i=1:length(lepesek)
    % Az 'integralo' programot lásd a dokumentum végén!
    hibak(i) = integralo(lepesek(i)); 
end
hibak

% Látható, hogy ahelyett, hogy egyre kisebb lépésközök esetén egyre kisebb hibákat kapnánk, egy idő után ezek növekedni fognak. Teljesül az alábbi állítás:
% Állítás: Egy másodrendű centrális közelítés esetén az optimális lépésköz
% h_{opt} = \sqrt[4]{\frac{48 \varepsilon}{M_4}},
% ahol \varepsilon a függvényértékek pontossága, és M_4 pedig a vizsgált függvény negyedik deriváltjának szuprémuma (a vizsgált intervallumon).
% Jelen esetben \sin^{(4)} (x) = \sin(x), azaz M_4 = \sup |\sin^{(4)}(x)| = \sup |\sin(x)| = 1, továbbá ismert, hogy a MATLAB esetén \varepsilon \approx 10^{-16}. Ezek alapján tehát
% h_{opt} = \sqrt[4]{\frac{48 \cdot 10^{-16}}{1}} \approx 2.6321 \cdot 10^{-4},
% tehát valóban 10^{-4} környékén lesz a legkisebb hiba.

h = nthroot(48*10^(-16),4)

integralo(h)

%
%% Numerikus integrálás
% Alapgondolat: Keressük az f(x) függvény határozott integrálját egy adott [a,b] intervallumon. Néha ismerjük f(x) alakját is, de gyakran csak az értékeit bizonyos pontokban.
% Híres közelítések:
% Trapézszabály:
% \int_a^b f(x) dx \approx |b-a| \left( \frac{1}{2} f(a) + \frac12 f(b)\right)
% Érintőformula:
% \int_a^b f(x) dx \approx |b-a| f\left( \frac{a+b}{2} \right)
% Simpson formula:
% \int_a^b f(x) dx \approx |b-a| \left( \frac16 f(a) + \frac46 f \left( \frac{a+b}{2}  \right)  + \frac16 f(b)\right)
% 3. feladat: Tekintsük az alábbi integrált:
% \int_0^2 \frac{1}{1+x^2} dx
% (a) Számítsuk ki az értékét a trapézmódszer és az érintőformula segítségével is! Mekkora lesz a hiba? (A pontos érték \arctan(2).)
% Megoldás: Jelen esetben f(x) = \frac{1}{1+x^2}. Alkalmazva a trapézszabályt:
% \int_0^2 \frac{1}{1+x^2} dx \approx 2 \left( \frac12 \frac{1}{1+0^2} + \frac12 \frac{1}{1+2^2}\right) = 1.2
% Hasonlóan az érintőszabály:
% \int_0^2 \frac{1}{1+x^2} dx \approx 2 \frac{1}{1+1^2} = 1.
% Mivel a pontos érték

atan(2)

% így a hibák:

abs(atan(2)-1.2)
abs(atan(2)-1)

% (b) Javítsuk az előző módszereket úgy, hogy ne a teljes [0,2] intervallumra alkalmazzuk őket, hanem számítsuk ki az integrálokat a [0,1] és [1,2] intervallumokon, majd az így kapott eredményket adjuk össze! (Ezeknek neve összetett trapéz, ill. összetett érintőformula.)
% Megoldás:  A trapézszabály:
% \int_0^2 \frac{1}{1+x^2} dx = \int_0^1 \frac{1}{1+x^2} dx + \int_1^2 \frac{1}{1+x^2} dx = \\
% \qquad  = 1 \left( \frac12 \frac{1}{1+0^2} + \frac12 \frac{1}{1+1^2}\right) + 1 \left( \frac12 \frac{1}{1+1^2} + \frac12 \frac{1}{1+2^2}\right) = \frac12 + \frac14 + \frac14 + \frac{1}{10} = 1.1
% Az érintőformula:
% \int_0^2 \frac{1}{1+x^2} dx = \int_0^1 \frac{1}{1+x^2} dx + \int_1^2 \frac{1}{1+x^2} dx = \\
% \qquad  = 1  \frac{1}{1+0.5^2} +  1 \frac{1}{1+1.5^2}  = \frac45 + \frac{4}{13} = \frac{72}{65} \approx 1.10769
% Most a hibák:

abs(atan(2)-1.1)
abs(atan(2)-72/65)

% (c) Vizsgáljuk meg, hogy mi történik, hogyha az imént a trapézszabályból származó x és az érintőformulából származó y értéket behelyettesítjük az (x+2y)/3 formulába! Mekkora lesz most a hiba? Alkalmazzuk a fenti módon a(z összetett) Simpson formulát is! Mit vehetünk észre?
% Megoldás:

uj = (1.1 + 2*72/65)/3;
abs(atan(2) - uj)

% A Simpson formulát használva a kisebb tartományokra:
% \int_0^2 \frac{1}{1+x^2} dx = \int_0^1 \frac{1}{1+x^2} dx + \int_1^2 \frac{1}{1+x^2} dx = \\
% \qquad  = 1 \left( \frac16 \frac{1}{1+0^2} + \frac46 \frac{1}{1+0.5^2} + \frac16 \frac{1}{1+1^2}\right) + 1 \left( \frac16 \frac{1}{1+1^2} + \frac46 \frac{1}{1+1.5^2} + \frac16 \frac{1}{1+2^2}\right) = \\
% \qquad =\frac16 + \frac{16}{30} + \frac{1}{12} + \frac{1}{12} + \frac{16}{78} + \frac{1}{30} \approx 1.105128
% Ekkor a hiba:

simpson = 1/6+16/30+2/12+16/78+1/30;
abs(atan(2) - simpson)

% A hiba ugyanakkora, hiszen a Simpson módszer pont a fenti módon származtatható a trapéz és érintő módszerekből.
%
% A 2. feladat megoldása:

function [hiba] = integralo(h)
kozelites = (sin(0.5-h) - 2*sin(0.5) + sin(0.5+h))/h^2;
hiba = abs(kozelites-(-0.479425538604203));
end