%% Hetedik gyakorlat
%% Nemlineáris egyenletek: Newton módszer.
%
% 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.
%
%% Egy változós eset
% Ötlet: Lépjünk a feltételezett zérushely irányába az érintő behúzásával, majd az így kapott ponton vizsgáljuk a függvény értékét.
%
% (Kép forrása: elméleti jegyzet)
% A módszer képlete:
% x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}.
% Tétel: Ha az f függvény "szép" (kétszer folytonosan deriválható), a vizsgált intervallumon első és második deriváltja se lesz sehol nulla és f''(x_0) f(x_0)>0, akkor a Newton-módszer konvergál a 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 az [1,2] intervallumon, í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
    f = x^3 + 2*x^2 + 10*x-20;
    x = x - 1/(3*x^2+4*x+10)*f
    
    k=k+1;
end

k

%
% Emlékeztető: Egy numerikus módszert p-ed rendűnek mondunk, ha az egymást követő közelítések e_k és e_{k+1} hibáira teljesül, hogy
% \lim_{k\to \infty \;} \;\frac{\|\;e_k \|\;}{\|\;e_{k+1} \|{\;}^p }=C
% valamely C nemnulla konstansra.
% 2. feladat: Keressük az alábbi függvény zérushelyét:
% f(x) = x^3 - 5 x^2 + 8 x - 4.
% Alkalmazzuk rá a Newton-iterációt! Honnan indítsuk, hogy konvergens legyen? Mit mondhatunk a módszer rendjéről?
% 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 konvergál az x=2 ponthoz, mivel a (2,5) intervallumon nem változik az előjel.
% A módszer rendje:

clear error
syms x  % szimbolikus számolásokhoz legyen x egy szimbolikus változó
f=x^3-5*x^2+8*x-4;
fder=diff(f); % kiszámítjuk az előző f függvény deriváltját (ehhez szimbolikus változóra van szükség)
f=matlabFunction(f); % átalakítjuk a szimbolikus függvényeket Matlab függvényekké 
fder=matlabFunction(fder);

x=5 
for i=1:15 
    i 
    x=x-f(x)/fder(x) 
    error(i)=abs(x-2);    
end
plot(log(error(1:end-1)),log(error(2:end)))
axis equal

% Látható, hogy a módszer csak elsőrendű (a vonal meredeksége a rend becslése). A gond, hogy a függvénynek x=2 kétszeres gyöke. Ahhoz, hogy a módszer ebben az esetben is másodrendű maradjon, használjuk a korábbiak helyett az x_{k+1} = x_k - m \frac{f(x_k)}{f'(x_k)} szabályt, ahol m a gyök multiplicitása. Ennek segítségével:

clear error

x=5 
for i=1:5 
    i 
    x=x-2*f(x)/fder(x)
    error(i)=abs(x-2);
end
plot(log(error(1:end-1)),log(error(2:end)))
axis equal

% Így a módszer már másodrendű.
%
%% Newton módszer alkalmazása egyenletrendszerekre
% Tekintsük az alábbi nemlineáris egyenletrendszert:
% f_1 (x_1, x_2, \dots x_n) = 0, \\
% f_2 (x_1, x_2, \dots x_n) = 0,\\
% \dots \\
% f_n (x_1, x_2, \dots x_n) = 0,
% ahol f_1,\; f_2,\dots \;f_n: \mathbb{R}^n \rightarrow \mathbb{R}  adott függvények.
% A Newton-módszer egyenletrendszerre:
% x_{k+1} = x_k - (J(x_k))^{-1} f(x_k),
% ahol (J(x))^{-1} az alábbi Jacobi mátrix inverze:
% J(x) = \left( \matrix{\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \dots & \frac{\partial f_1}{\partial x_n} \cr
% \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2 }{\partial x_2} & \dots & \frac{\partial f_2}{\partial x_n} \cr \vdots & \vdots & \vdots & \vdots \cr \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \dots & \frac{\partial f_n}{\partial x_n}} \right).
% Megfelelő feltételek mellett a módszer konvergens lesz.
%
% 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