%% Hatodik gyakorlat
%% Nemlineáris egyenletek: intervallumfelezés, fixpont iteráció
%
% Cél: Az f(x)=0 egyenlet megoldása, ahol f adott, szép tulajdonságú (pl. folytonosan deriválható) függvény, azaz keressük az f függvény gyökeit/zérushelyeit. A mai gyakorlaton kizárólag f(x): \mathbb{R} \rightarrow \mathbb{R} függvényeket vizsgálunk.
%
%% Intervallumfelezés
% Tegyük fel, hogy tudjuk, hogy egy adott [a,b] intervallumon a függvényünknek van egy gyöke: például ezt onnan tudhatjuk, hogy f(a)<0 és f(b)>0. Az általános eset megszorítása nélkül ez feltehető (a másik esetben hasonlóan megy az algoritmus).
% Alapgondolat: Vizsgáljuk meg a függvény előjelét az \frac{a+b}{2} pontban! Ha itt pozitív, az azt jelenti, hogy a keresett zérushely az \left[a,\; \frac{a+b}{2}\right) intervallumban van, azaz itt érdemes folytatni a keresést. Ugyanakkor ha itt negatív, akkor a zérushely az \left( \frac{a+b}{2},\; b\right] intervallumban található, és itt kell tovább keresnünk. Ezt követően az előbbi gondolatot folytatjuk az új, kisebb intervallumokon, azaz vagy az \frac{a+b}{4}, vagy a \frac{3 (a+b)}{4} pontban vizsgáljuk a előjelet, és így tovább.
% Az algoritmus addig fut, amíg f értéke a vizsgált pontban nulla lesz (ez programok esetén a gépi hibák miatt nem reális), vagy amikor a vizsgált intervallum már nagyon kicsi ("kvázi nulla hosszúságú"). Természetesen megadhatunk maximális lépésszámot is, amely után álljon le mindenképpen a program.
% % tol: az elvárt pontosság.
%
% while (b - a)>tol
%      x = (a+b)/2;
%      f = f(x);
%
%      if f == 0
%         return x; % Ha véletlenül pont a felezőponton van a zérushely, akkor megállunk
%      else
%         if f>0
%              b = x;
%         else
%              a = x;
%         end
%      end
% end
% A módszer mindenképpen megtalál egy zérushelyet (ha több is van, akkor a zérushelyek közül valamelyiket).
%
% 1. feladat: Fibonacci egyik cikkében szerepel az alábbi állítás: az
% f(x) = x^3 + 2 x^2 + 10 x - 20
% egyenlet gyöke x = 1.368808107... Ellenőrizzük le az intervallumfelezéses módszerrel, hogy jól számolt-e!
% Megoldás: Legyen a toleranciánk 10^{-10}, és indítsuk az iterációt az [1,2] intervallumon!

format long
tol=10^(-10);

a=1;
b=2;

k=0;

while abs(b-a)>tol
    x = (a+b)/2
    f = x^3+ 2*x^2+ 10*x-20;
    if f==0
        return 
    else
        if f>0
            b=x;
        else
            a=x;
        end
    end
    k=k+1;
end
x
k

% Látható, hogy 34 lépés kellett a mószernek, de megtalálta a zérushelyet.
% A függvényt is ábrázolhatjuk:

x = linspace(-1,5,1000);
f = x.^3+ 2*x.^2+ 10*x-20;
plot(x,f)

%
%% Fixpont iteráció
% Alapgondolat: A korábbi f(x)=0 egyenletünket átírhatjuk
% x + c f(x) = x
% alakba, ahol c egy, a módszer közben meghatározandó konstans. Vegyük észre, hogy így a feladat tulajdonképpen az F(x) = x+ c f(x) függvény fixpontjának megkeresésével ekvivalens.
% Másik módszer: Ha az f(x) függvényből kifejezhető x, azaz pl. f(x) = k(x)  + x\;g(x), akkor egy másik lehetőség egy fixpont-egyenlet felírására az alábbi:
% x = - \frac{k(x)}{g(x)}
% Ekkor tehát az F(x) = - \frac{k(x)}{g(x)} függvény fixpontját keressük.
% Tétel: Ha F függvény ráképezés (tehát ha  [a,b]-n van értelmezve, akkor értékei is az [a,b] intervallumban vannak) és kontrakció a vizsgált [a,b] intervallumon, akkor az x_{k+1}=F(x_k) iteráció tart az F függvény fixpontjához.
% Állítás: Ha F függvény folytonos [a,b]-n, differenciálható (a,b)-n és van egy olyan 0 \leq q < 1 érték, amelyre |F'(x)|< q, akkor F kontrakció.
% 2. feladat: Keressük meg az első feladatban szereplő
% f(x) = x^3 + 2 x^2 + 10 x - 20
% egyenlet gyökét fixpont iteráció segítségével is! Legyen a vizsgált intervallum [a,b]=[1,2], és addig lépjünk, amíg a hiba kisebb, mint 10^{-10}.
% Megoldás (első módszer): Tekintsük az alábbi egyenletet:
% x = x + c \;(x^3 + 2 x^2  + 10x - 20).
% Tehát legyen F(x) =  x + c \;(x^3 + 2 x^2  + 10x - 20). Válasszuk úgy a c értéket, hogy teljesüljenek a tétel feltételei, azaz hogy F kontrakció és ráképezés legyen. Az elsőhöz számítsuk ki a deriváltat:
% |F'(x)| = |1 + 3 c x^2 + 4 c x + 10 c|< 1
% Látható, hogy c=-0.01 esetén a függvény kisebb, mint 1.

x = linspace(0,3);
c= -0.01;
f = abs(1+3*c*x.^2+ 4*c*x+ 10*c);
plot(x,f)

% Ennek precíz bizonyítása: vizsgáljuk meg az 1 + 3 c x^2 + 4 c x + 10 c függvényt! Ennek szélsőértéke a 6 cx + 4c=0 egyenlet megoldása, ami c \neq 0 esetén x = \frac46, ami kívül esik a vizsgált [1,2] intervallumon, tehát a deriváltfüggvény monoton itt, így értékeinek vizsgálatához elegendő a két végpontban megvizsgálni az értékét: az x=1-ben ez 1 + 3c + 4c + 10 c = 1 + 17 c, míg az x=2-ben 1 + 12 c + 8 c + 10 c = 1 + 30 c. Azt szeretnénk, ha mindkettő kisebb lenne abszolútértékben, mint 1. Látható, hogy csak a c<0 esetben lehetnek kisebbek, mint 1, és c>-\frac{2}{30} \approx -0.067 esetén lesznek nagyobbak, mint -1. Így tehát a c=-0.01 jó választás.
% A ráképezéshez arra van szükség, hogy 1 \leq F(x) \leq 2 igaz legyen minden x \in [1,2] értékre. Látható, hogy a c=-0.01 itt is jó választás.

x = linspace(0,3);
c= -0.01;
f = x + c*(x.^3 + 2*x.^2+ 10*x-20);
plot(x,f)

% Ennek precíz bizonyítása: vizsgáljuk meg az x + c \;(x^3 + 2 x^2  + 10x - 20) függvény monotonitását az [1,2] intervallumon! Ennek deriváltja (az előtő lépésben kiszámoltuk): 1 + 3 c x^2 + 4 c x + 10 c, így a szélsőértékek -\frac23 \pm \frac{\sqrt{-26 c^2  - 3 c}}{3c}. Ha használjuk az előző lépésben használt c=-0.01választást, akkor a gyökök kívül esnek a [0,1] intervallumon, így itt biztosan monoton a függvény. Ezért elegendő az értékeit kiszámítani a határpontokban: F(1)= 1 + c + 2 c + 10 c - 20 c = 1 - 7 c, ami a fenti választás mellett nagyobb, mint 1. Hasonlóan, F(2) = 2 + 8 c + 8 c + 20 c - 20 c = 2 + 16c, ami a fenti választás mellett kisebb, mint 2. Ezek alapján tehát valóban ráképezést kaptunk.
% Mindezek alapján tehát a c = -0.01 választás mellett konvergens lesz a módszer.

x = 1.3;
x_elozo=1000;
c=-0.01;
k=0;

while abs(x-x_elozo)>10^(-10)
    x_elozo=x;
    x = x + c*(x^3 + 2*x^2+ 10*x-20) 
    k=k+1;
end
k

% Megoldás (második módszer): Rendezzük át az eredeti egyenletet az alábbi módon:
% x = \frac{20}{x^2+2x+10}.
% Tehát legyen F(x) = \frac{20}{x^2+2x+10}. Erről kell belátni, hogy kontrakció és ráképezés. Ezen függvény és deriváltja is monoton az [1,2] intervallumon. Az első feltételhez:
% |F'(x)| = \left| \frac{-20(2x+2)}{(x^2+2x+10)^2} \right| \leq 20 \frac{6}{13^2} = \frac{120}{169}<1,
% ahol felhasználtuk, hogy x \in [1,2]. Hasonlóan, a ráképezésnél:
% 1 \leq \frac{20}{18} \leq \frac{20}{x^2+2x+10}\leq \frac{20}{13}\leq 2,
% tehát valóban ráképezés.

x = 1;
x_elozo=1000;
k=0;

while abs(x-x_elozo)>10^(-10)
    x_elozo=x;
    x = 20/(x^2+2*x+10)
    k=k+1;
end
k

%
% 3. feladat: Határozzuk meg azon pontot, ahol az f(x)=x^5 - x függvény értéke 2 lesz. Használjunk intervallumfelezést és fixpont iterációt is!
% Megoldás (intervallumfelezés): Azaz most az f(x)=x^5 - x - 2 függvény zérushelyét keressük. Behelyettesítve néhány értéket kapjuk, hogy f(1) = -2 és f(2) = 28, azaz érdemes az [1,2] intervallumon keresni a gyököt.

format long
tol=10^(-10);

a=1;
b=2;

k=0;

while abs(b-a)>tol
    x = (a+b)/2
    f = x^5-x-2;
    if f==0
        return 
    else
        if f>0
            b=x;
        else
            a=x;
        end
    end
    k=k+1;
end
x
k

% Megoldás (fixpont): Azaz most az f(x)=x^5 - x - 2 függvény zérushelyét keressük. Ehhez fel kellene írnunk egy fixpont iterációt. Behelyettesítve néhány értéket kapjuk, hogy f(1) = -2 és f(2) = 28, azaz érdemes az [1,2] intervallumon keresni a gyököt.
% Első ötlet: Átrendezve x = x^5 - 2 adódik, azaz legyen F(x) = x^5 - 2. Ekkor azonban F'(x) = 5 x^4 , ami nem lesz kisebb, mint 1 az [1,2] intervallumon, így nem lesz jó.
% Második ötlet: Fejezzük ki az egyenletben szereplő másik x-et! Ekkor x = \sqrt[5]{x + 2} adódik, azaz legyen F(x) = \sqrt[5]{x + 2}. Ekkor
% |F'(x)| = \left| \frac15 \left( x+2 \right)^{-4/5} \right| = \left| \frac15 \frac{1}{(x+2)^{4/5}} \right| \leq 0.09 < 1,
% azaz kontrakció lesz (használtuk, hogy a függvény monoton az [1,2] intervallumon). A ráképezés az alábbi becslésből adódik:
% 1 \leq \sqrt[5] 3 \leq \sqrt[5]{x + 2} \leq \sqrt[5] 4 \leq 2.
% Azaz így mivel ez egy ráképezés és kontrakció, így a fixpont iteráció konvergálni fog.

format long
x = 1;
x_elozo=1000;
k=0;

while abs(x-x_elozo)>10^(-10)
    x_elozo=x;
    x = nthroot(x+2,5)
    k=k+1;
end
k