%% Nyolcadik gyakorlat
%% Nemlineáris egyenletek: Newton-módszer (Newton-iteráció)
%
% Cél: Az f(x)=0 egyenlet(rendszer) 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.
% Ötlet 1: Ha x egy pont, x^* pedig az f függvény zérushelye, akkor
% f(x) = f(x) - f(x^*) \approx f'(x)(x - x^*)
% és innen x^* egy közelítése kifejezhető:
% x^* \approx x - [f'(x)]^{-1} f(x),
% így legyen a módszer képlete
% x_{k+1} = x_k - [f'(x_k)]^{-1} f(x_k).
% egyváltozós esetben ezt úgy is megfogalmazhatjuk, hogy az x_k pontba behúzzuk a függvény érintőjét és megkeressük, hogy ez az egyenes hol metszi az x tengelyt; legyen ez a pont az x_{k+1}.
%
% (Kép forrása: elméleti jegyzet)
% Ötlet 2: Már korábban is láttuk, hogy az f(x) = 0 egyenletből képzett x \mapsto x + f(x) függvény általában nem lesz kontrakció, hiába lenne ennek fixpontja az előbbi egyenlet megoldása. Ehelyett valamilyen egyéb átalakításokkal szoktunk próbálkozni.
% Korábban szerepelt az f(x) = 0 \iff cf(x) = 0 ötlet, ahol c egy nem nulla konstans. Ha f lineáris; vagy nemlineáris, de monoton, akkor láttuk, hogy ez a megközelítés nem reménytelen. Ebben a két esetben azt is láttuk, hogy a c konstans megválasztásánál fontos szerepet kaptak f' (pontonkénti) sajátértékei.
% Bizonyos értelemben a másik véglet az
% f(x) = 0 \iff f^{-1}(f(x)) = f^{-1}(0) \iff x = f^{-1}(0).
% ötlet, amivel az egyenlet baloldalán könnyen kezelhető függvényt kapunk, de amennyivel a baloldal egyszerűsödik, a jobboldal annyival bonyolódik.
% Ezen két véglet között helyezkedik el az
% f(x) = 0 \iff [f'(x)]^{-1}\cdot f(x) = [f'(x)]^{-1} \cdot 0 = 0
% átalakítás, mely ekvivalens, ha f'(x) invertálható (a nulla nem sajátértéke). Ezzel az ötlettel kapjuk az
% F(x) = x - [f'(x)]^{-1}f(x)
% függvényre alkalmazható fixpont-iterációt, vagyis a Newton-iterációt.
% Megjegyzés: Az iteráció átírható az
% \begin{cases}
% f'(x_k) \, d_k &= - f(x_k),\\
%  x_{k+1}  &= x_k + d_k.
% \end{cases}
% alakba is, amiből talán még egyértelműbb, hogy ezzel a módszerrel egy nemlineáris egyenlet(rendszer) megoldását vezetjük vissza (f'(x_k) együtthatómátrixú) lineáris egyenlet(rendszer)ek egymásutáni megoldására.
% Ezen a ponton azt is meggondolhatjuk, hogy ha a lineáris esetben -azaz haf(x) = Ax - b- a módszer egy lépéssel megtalálja az Ax = b egyenlet megoldását.
% Rend: A módszer sajátossága, hogy azt a zérushely kellően kicsi környezetéből indítva viszonylag természetes feltételek mellett garantálható, hogy a konvergencia másodrendű lesz. Ennek megértéséhez vizsgáljuk meg a módszer által generált hibavektorokat! Legyen tehát
% x_{k+1} = x_k - [f'(x_k)]^{-1} f(x_k),
% ekkor x_{k+1} -x^*  =
%                          = x_k-x^* - [f'(x_k)]^{-1} \left(f(x_k) - f(x^*)\right) \\
% = [f'(x_k)]^{-1} \left( f'(x_k)(x_k-x^*) + f(x^*) - f(x_k)\right) \\
% = [f'(x_k)]^{-1} \left(- f'(x_k)(x^*-x_k) + \int_0^1f'(x_k + t(x^*-x_k))(x^*-x_k)\, dt\right) \\
% = [f'(x_k)]^{-1} \int_0^1\left( f'(x_k + t(x^*-x_k)) -  f'(x_k) \right)(x^*-x_k)\, dt.
% Ezért ha x^* egy gömbi környezetén tetszőleges x,y pontokra teljesül, hogy
% \| [f'(x)]^{-1}\| \leq \frac {1}{h} és
% \|f'(x) - f'(y)\| \leq L \|x-y\|,
% akkor a fentiekből kapjuk, hogy
% \|x_{k+1} -x^*\| \leq \| [f'(x_k)]^{-1}\|  \|x_k - x^*\|^2 L\int_0^1 t\, dt \leq \frac{L}{2h} \| x_k - x^* \|^2,
% amiből láthatjuk, hogy a módszerrel az x^* zérushelyhez kellően közelről indulva ezen feltételek mellett a hibák másodrendben csökkennek.
%
%% Egyváltozós esetben,
% ha a gyököt egy adott [a,b] intervallumban keressük, akkor kimondhatjuk a következő tételt a módszer konvergenciájával kapcsolatban:
% Tétel: Ha egy [a,b] intervallumon teljesül, hogy azon
% f kétszer folytonosan differenciálható,
% se az f', se az f'' függvénynek nincs gyöke, és
% egy x_0 pontban f''(x_0) f(x_0)>0,
% akkor az x_0 pontból indított Newton-módszer konvergál az x^* zérushelyhez. Ha pedig ezeken túl az is teljesül, hogy az intervallumon
%  m_1\leq |f'| \leq M_1 \qquad \text{és} \qquad  m_2 \leq |f''| \leq M_2,
% valamely 0 < m_i \leq M_i konstansokkal; továbbá az x_0 kezdőpontra
% |x_0 - x^*| < \mathrm{min} \left\{ |a-x^*|,\; |b-x^*|, \;\frac{2 m_1}{M_2} \right\},
% akkor a módszer másodrendben konvergál az x^* zérushelyhez.
%
% 1. feladat: Keressük meg az előző gyakorlaton bemutatott
% f(x) = x^3 + 2 x^2 + 10 x - 20
% függvény zérushelyét a Newton-módszer segítségével!
% Megoldás: Az előző gyakorlaton már láttuk, hogy f(1) <0 és f(2)>0, azaz érdemes az [a,b]=[1,2] intervallumon vizsgálódni. Az indítási hely megállapításához számítsuk ki az f
% függvény deriváltjait!
% f'(x) = 3 x^2 + 4 x + 10
% f''(x) = 6 x + 4
% Ez mindenhol pozitív, így a f''(x_0) f(x_0)>0 feltétel teljesüléséhez pl. az x_0=2 kezdőpont megfelelő lesz. Mivel mindkét derivált pozitív, így a tétel további feltételei is teljesülnek, azaz konvergens lesz a módszer.

format long

x = 2;
tol = 10^(-10);
f = x^3 + 2*x^2 + 10*x-20;
k=0;

while f>tol
    x = x - 1/(3*x^2+4*x+10)*(x^3 + 2*x^2 + 10*x-20)
    f = x^3 + 2*x^2 + 10*x-20;
    k=k+1;
end

k

%
% 2. feladat: Keressük az alábbi függvény zérushelyét:
% f(x) = x^3 - 5 x^2 + 8 x - 4.
% (a) Alkalmazzuk rá a Newton-iterációt! Honnan indítsuk, hogy konvergens legyen?
% Megoldás: Látható, hogy f(0)=-4 és f(5)=36>0, azaz itt biztosan lesz egy zérushelye. Vizsgáljuk meg a deriváltakat x=5-ben!
% f'(x) = 3 x^2-10x + 8
% f''(x) = 6 x - 10,
% azaz ezek alapján f(5) f''(5) = 36 \cdot 20>0. Látható, hogy az első derivált negatívvá válik x \in (4/3, 2) esetén, továbbá a második derivált is negatív lesz x<4/3 esetén. Ennek ellenére nézzük meg, mi történik a módszerben:

format long

x = 5;
tol = 10^(-10);
f = x^3 - 5*x^2 + 8*x-4;
k=0;


while abs(f)>tol
    x = x - 1/(3*x^2-10*x+8)*(x^3 - 5*x^2 + 8*x-4)
    f = x^3 - 5*x^2 + 8*x-4;
    k=k+1;
end
k

% Tehát a módszer a feltételnek megfelelően konvergál az x=2 ponthoz, mivel a (2,5) intervallumon nem változik az előjel.
% (b) Vizsgáljuk meg a módszer rendjét! Mit tapasztalunk? Indokoljuk a jelenséget!
% Megoldás: Az előző feladatsoron szereplő definíciót felhasználva:

format long

x = 5;
tol = 10^(-10);
f = x^3 - 5*x^2 + 8*x-4;
rendek=zeros(20,1);

for i=1:20
    x_elozo=x;
    x = x - 1/(3*x^2-10*x+8)*(x^3 - 5*x^2 + 8*x-4);
    f = x^3 - 5*x^2 + 8*x-4;
    hiba=abs(x-2);
    hiba_elozo = abs(x_elozo-2);
    rendek(i)=log(hiba)/log(hiba_elozo);
end

rendek

% A módszer csak elsőrendű. Megvizsgálva az előző tétel feltételeit láthatjuk, hogy
% |a-x^*| = |0-2|=2,
% |b - x^*|=|5-2|=3,
% m_1 = \mathrm{min} |f'(x)| = 0,
% M_2 = \max |f''(x)| = 20.
% Mindezek alapján tehát |x_0 - x^*| \leq 0 feltételnek kellene teljesülnie, ami persze nem lesz igaz, így a tétel nem garantálhat nekünk másodrendű konvergenciát. A probléma az, hogy |f'(x)| értéke x \rightarrow x^* esetén nullához konvergál, azaz nincs pozitív alsó korlátja.
%
%% Többváltozós esetben
% egy-egy egyenletrendszer megoldásáról (ekvivalensen: \mathbb{R}^n-beli gyökkeresési feladatról) beszélhetünk. A rendszert alkotó egyenleteket az
% f_1 (y_1, y_2, \dots y_n) = 0, \\
% f_2 (y_1, y_2, \dots y_n) = 0,\\
% \dots \\
% f_n (y_1, y_2, \dots y_n) = 0,
% nullára rendezett alakba írjuk, ahol f_1,\; f_2,\dots \;f_n: \mathbb{R}^n \rightarrow \mathbb{R}  adott függvények, melyek pontosan az f:\mathbb{R}^n \to \mathbb{R}^n függvény komponensei lesznek.
% A Newton-módszerben szereplő [f'(y)]^{-1}  kifejezés ekkor nem más, mint az alábbi Jacobi-mátrix inverze:
% f'(y) = J(y) = \left( \matrix{\frac{\partial f_1}{\partial y_1} & \frac{\partial f_1}{\partial y_2} & \dots & \frac{\partial f_1}{\partial y_n} \cr
% \frac{\partial f_2}{\partial y_1} & \frac{\partial f_2 }{\partial y_2} & \dots & \frac{\partial f_2}{\partial y_n} \cr \vdots & \vdots & \vdots & \vdots \cr \frac{\partial f_n}{\partial y_1} & \frac{\partial f_n}{\partial y_2} & \dots & \frac{\partial f_n}{\partial y_n}} \right).
%
% 3. feladat: Tekinsük az alábbi egyenletrendszert:
% \begin{cases}
% 5 x^2 - y^2 &= 0,\\
% y - 0.25(\sin x +\cos y) &= 0.
% \end{cases}
% Newton-módszer segítségével keressük ennek a (x_0,y_0)=(0.25,0.25) pont környezetében lévő megoldását!
% Megoldás: Ki kell számítanunk a Jacobi mátrixot (ahol f_1(x,y) = 5x^2 - y^2 és f_2(x,y) = y - 0.25(\sin x +\cos y) = 0):
% J(x) = \left( \matrix{ 10 x & -2y \cr
% -0.25 \cos x & 1 + 0.25 \sin y } \right).
% Így a kód:

format long

x = [0.25 0.25]';
tol = 10^(-10);
f = [5*x(1)^2-x(2)^2  x(2)-0.25*(sin(x(1))+cos(x(2)))];
k=0;


while norm(f,'Inf')>tol
    J = [10*x(1) -2*x(2); -0.25*cos(x(1)) 1+0.25*sin(x(2))];
    f = [5*x(1)^2-x(2)^2; x(2)-0.25*(sin(x(1))+cos(x(2)))];
    x = x - inv(J)*f
    k=k+1;
end
k

% Használhatjuk a lenti kódot is:

f = @(x, y)([                   ...
    5*x^2-y^2;                  ...
    y-0.25*(sin(x)+cos(y))      ...
]);
df = @(x, y)([                  ...
         10*x,          -2*y;   ...
 -0.25*cos(x), 1+0.25*sin(y)    ...
]);

x0 = [1; 1]/4;

[r, num_steps] = newton_it(     ...
    @(x)( f(x(1), x(2))),       ...
    @(x)(df(x(1), x(2))),       ...
    x0,                         ...
    1e-10, 100                  ...
)

%
% Megjegyzés: Inegzakt- vagy Kvázi-Newton módszerek
% Láthattuk, hogy a Newton-iteráció viszonylag természetes feltételek mellett másodrendben konvergál a zérushelytől nem túl messziről indítva, azonban abból is kaphattunk egy ízelítőt, hogy ennek ára van: lépésenként ki kell értékelni a függvény deriváltját és ezzel egy-egy lineáris egyenletet is meg kell oldani. Ez főleg nagyobb rendszerek esetén praktikusan reménytelenül számításigényes feladat. Egy lehetséges feloldása ennek a problémának, ha az
% f'(x_k) \, d_k = - f(x_k)
% egyenletet nem mindig oldjuk meg pontosan és d_k helyett egy \tilde d_k vektorral lépünk arrébb. Pontosabban:
% \begin{cases}
% \| f'(x_k) \, \tilde d_k + f(x_k)\| & \text{legyen kellően kicsi},\\
%  x_{k+1}  &= x_k + \tilde d_k.
% \end{cases}
% Egy egyszerű alesete ennek, ha azf'(x_k) kifejezést csak néhány lépésenként értékeljük ki, egyébként mindig a legutóbbi mátrixszal oldunk meg egyenletet. Mivel feltesszük, hogy az f' függvény Lipschitz-folytonos, ezért ha nem lépünk nagyot (azaz \|x_k - x_{k+1}\|kicsi), akkor intuitíve ezzel olyan nagyot nem hibázunk. Ha az így kapott mátrix LU-felbontását kiszámoljuk, majd ezt használjuk az egyenletmegoldáshoz, akkor ezzel is időt spórolunk.
%
% Szorgalmi: Oldjuk meg a 3. feladatot úgy, hogy csak k lépésenként értékeljük ki (újra) a J(x) mátrixszot. Kiértékelés után készítsük is el egyből a mátrix LU-felbontását, és használjuk az így kapott mátrixokat a következő k lineáris egyenletrendszer megoldására. Különböző k értékek mellett (pl. k=1,2,3,4,5) becsüljük meg a módszer rendjét és mérjünk futásidőt.

%

%
% Szorgalmi: A 2. feladatban láttuk, hogy ha az
% f(x) = x^3 - 5 x^2 + 8 x - 4
% függvény zérushelyét keressük Newton módszer segítségével az x_0=5 pontból indulva, akkor a módszer csupán elsőrendű lesz.
% (a) Minden Newton módszer felfogható úgy, mint egy fixpont iteráció az F(x) = x - \frac{f(x)}{f'(x)} választással. Lássuk be az előző gyakorlaton tanult tétel segítségével, hogy a jelen feladatban szereplő iteráció valóban nem lehet másod, csak elsőrendű! (Tipp: ha nem értelmes F'(x^*), vizsgáljuk \lim_{x \rightarrow x^*} F'(x) értékét!)

%

% (b) Az előző gondolatmenetet folytatva módosítsuk a Newton iterációt úgy (a \frac{f(x)}{f'(x)} tag megfelelő konstanssal való szorzásával), hogy a módszer már másodrendű legyen! Ellenőrizzük is ezt le numerikusan!

%

%
% Függvénytemető:

[r, num_steps] = newton_it(@(x)(x*exp(-x^2/2)), @(x)((1-x^2)*exp(-x^2/2)), 1/4, 1e-10, 100)
[r, num_steps] = newton_it(@(x)(x), @(x)(eye(numel(x))), ones(10,1), 1e-10, 100)

%

function [x, num_steps] = newton_it(f, df, x0, atol, max_steps)
    x = x0;
    for num_steps = 1:max_steps
        d = - df(x) \ f(x);
        x = x + d;
        if norm(d) < atol
            break;
        end
    end
end