%% Tizedik gyakorlat
%% Köz. diff. egyenlet numerikus megoldása: Euler módszer
%
% Kedvcsináló: https://www.youtube.com/watch?v=v-pbGAts_Fg
%
% Cél: Szeretnénk megoldani az alábbi közönséges differenciálegyenletet:
% \begin{cases}
% u'(t) = f(t, u(t)), \\
% u(0) = u_0, \qquad \quad \;
% \end{cases}
% ahol f olyan "szép" függvény, hogy ennek létezik egyértelmű megoldása, valamint u_0 adott kezdeti érték. Ha f nem speciális alakú, akkor nem ismert általános eljárás, ami zárt alakban visszaadná a megoldást, ezért numerikusan közelítjük a megoldást. Gyakran nem is kell maga a függvény pontos formulája, csak a függvény értéke valamilyen t = T esetén.
%
% Euler módszer: Közelítsük a keresett u(T) értéket úgy, hogy felbontjuk a [0,T] tartományt kis intervallumokra - legyen ezen intervallumok hossza h, és az így kapott pontok pedig t_0=0, \;\; t_1=h, \;\; t_2 = 2h \;\; \dots \;\; t_N = T. Legyen egy adott pontban a közelítés u(t_n)\approx y_n. Belátható (az 5. gyakorlaton tanult módszerekkel), hogy
% u'(t_n)\approx \frac{u(t_{n+1}) - u(t_n)}{h}
% és
% u'(t_n) \approx \frac{u(t_n) - u(t_{n-1})}{h}.
% A fentieket beírva a korábbi egyenletbe kapjuk az alábbi két algorimtust:
% \frac{y_{n+1} - y_{n}}{h} = f(t_n, y_n)
% \frac{y_{n} - y_{n-1}}{h} = f(t_n, y_{n})
% Ha ezeket átrendezzük, kapjuk az explicit Euler módszert:
% y_{n+1} = y_n + h f(t_n, y_{n}),
% illetve az implicit Euler módszert:
% y_n = y_{n-1} + h f(t_n, y_n).
% Látható, hogy az explicit esetében könnyedén kiszámítható mindig a következő elem, míg az implicit esetén minden lépésben meg kell oldanunk ezt az egyenletet valamilyen módszerrel (pl. Newton módszerrel).
% 1. feladat: Oldjuk meg az alábbi egyenletet explicit Euler módszer segítségével a [0,10] intervallumon! Hogyan néz ki ennek egy lépése az implicit módszerrel? Oldjuk meg így is Newton iterációt alkalmazva (ez addig fusson, amíg az egymás követő értékek különdsége már kisebb, mint 10^{-8}, és a kezdeti érték legyen mindig az előző lépésben kiszámolt y)!
% \begin{cases}
% y'(t) = \mathrm{arctg} (y(t)), \\
% y(0) = 1. \qquad\quad\;\;\;\;
% \end{cases}
% Megoldás: Ez egy szétválasztható egyenlet, de a megoldása során az \displaystyle\int \frac{1}{\mathrm{arctg}(y)} dy integrált kellene kiszámolni, ami nehéz (nincs is zárt alakja a függvénynek).
% Az explicit Euler módszer egy lépése:
% y_n = y_{n-1} + h \;\mathrm{arctg} (y_{n-1})
% Kódját lásd a dokumentum végén, ábrázolás h=5 és h=0.01 esetén:

[durvaido, durvaertek]=elso(5);
[pontosido, pontosertek]=elso(0.01);

plot(durvaido, durvaertek, pontosido, pontosertek)
legend('durva','pontosabb')

% Az implicit Euler módszer képletét ha felírnánk:
% y_n = y_{n-1} + h \; \mathrm{arctg(y_n)}
% Átrendezve:
% y_n - h\; \mathrm{arctg}(y_n) = y_{n-1}
% Látható, hogy ehhez az g(x)=x - h\; \mathrm{arctg} (x) függvényt kellene invertálni. Ehelyett valójában az alábbi függvény zérushelyét keressük minden lépésben:
%  f(x)=x - h\; \mathrm{arctg}(x) - y_{n-1}.
% Így a konkrét Newton iteráció alakja:
% x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} = x_k - \frac{x_k - h \text{arctg}(x_k) - y_{n-1}}{1-\frac{h}{x_k^2+1}}.

% Lépésköz
h=0.5;

% Végső érték
T=10;

% Lépések száma a [0,T] intervallumon (a kezdetit is beleszámolva)
N = T/h+1;

% Felveszünk egy vektort, amiben tároljuk a lépéseket.
ertekek = zeros(N,1);

% Kezdeti feltétel:
ertekek(1)=1;

% Egy másikban tároljuk az időket (az ábrázoláshoz kell ez majd)
idok = zeros(N,1);

for k=2:N
    x_elozo=10^6;
    x=ertekek(k-1);
    while abs(x-x_elozo)>10^(-2)
        x_elozo=x;
        x=x-(x-h*atan(x)-ertekek(k-1))/(1-h/(x^2+1));  
    end
    ertekek(k) = x;
    idok(k)=h*(k-1);
end

plot(idok,ertekek)

%
% 2. feladat: Tekintsük az alábbi egyenletet:
% \begin{cases}
% u'(t) = - 100 \;u(t), \\
% u(0) = 1. \qquad\quad\;\;
% \end{cases}
% Számítsuk ki a pontos megoldást, majd vizsgáljuk meg, hogy mi történik, hogyha az explicit, illetve, ha az implicit módszert alkalmazzuk az egyenletre a [0,1] intervallumon h = 0.2 lépésközzel! Indokoljuk meg a kapott eredményt!
% Megoldás: A pontos megoldás u(t) = e^{-100 t}, ennek értéke u(1) = e^{-100} \approx 3.72 \cdot 10^{-44}.
% Felírva az explict módszert:
% y_n = y_{n-1} + h (-100) y_{n-1} = (1-100h) y_{n-1} = (1-100h)^n \; y_0 = (1-20)^n \; y_0 = (-19)^n \; y_0
% MATLAB-ban ezeket kiszámítva:


% Időértékek
t = 0;
h=0.2;
y=zeros(1/h+1,1);

% Kezdeti érték
y(1) = 1;


i=2;
while t<1
    y(i) = y(i-1) - 100*h*y(i-1);
    t = t + h;
    i=i+1;
end

plot(y)

% Látható, hogy nagyon rosszul viselkedik a módszer (nem stabil, lásd előadás).
% Az implicit Euler módszert felírva:
% y_n = y_{n-1} + h (-100) y_n.
% Átrendezve:
% (1+100h) y_n = y_{n-1},
% y_{n} = \frac{1}{1+100h} y_{n-1} = \left(\frac{1}{1+100h}\right)^n y_0 = \frac{1}{21^n} y_0.
% Ez már szépen viselkedik:

% Időértékek
t = 0;
h=0.2;
y=zeros(1/h+1,1);

% Kezdeti érték
y(1) = 1;

i=2;
while t<1
    y(i) = 1/(1+100*h)*y(i-1);
    t = t + h;
    i=i+1;
end

plot(y)

% Megmutatható, hogy y_n = \left( \frac{1}{\frac{100}{n+1} + 1}\right)^{n+1} \rightarrow e^{-100} ha n \rightarrow \infty, azaz egyre kisebb lépésközök mellett a módszer valóban konvergál (stabil).
%
% 3. feladat: Oldjuk meg explicit Euler módszerrel az alábbi egyenletrendszert a t \in [0,1] intervallumon h=0.2-vel és u_1(0)=2, u_2 (0)=-3 kezdeti értékekkel!
% \begin{cases}
% u_1'(t) = 8 u_1(t) + 9 u_2(t),\;\;\;\;\;\; \\
% u_2'(t) = -18 u_1(t) - 19 u_2(t).
% \end{cases}
% Megoldás: Az Euler módszert alkalmazzuk egy rendszerre: ha az egyenletek jobb oldalait f_1(t, u_1,u_2)-vel és f_2(t,u_1, u_2)-vel jelöljük, akkor a módszer:
% \left( \begin{array}{c}
% y_{1,n}\\
% y_{2,n}
% \end{array}\right) = \left( \begin{array}{c}
% y_{1,n-1}\\
% y_{2,n-1}
% \end{array}\right)  + h \left( \begin{array}{c}
% f_1(t_{n-1},y_{1,n-1}, y_{2,n-1})\\
% f_2(t_{n-1},y_{1,n-1},y_{2,n-1})
% \end{array}\right).
% Leprogramozva:

% Lépésköz
h=0.001;

% Végső érték
T=1;

% Lépések száma a [0,T] intervallumon (a kezdetit is beleszámolva)
N = T/h+1;

% Felveszünk egy vektort, amiben tároljuk a lépéseket.
ertekek = zeros(N,2);

% Kezdeti feltétel:
ertekek(1,1)=2;
ertekek(1,2)=-3;

% Egy másikban tároljuk az időket (az ábrázoláshoz kell ez majd)
idok = zeros(N,1);

for k=2:N
    ertekek(k,:) = ertekek(k-1,:) + h*([8 9; -18 -19]*ertekek(k-1,:)')';
    idok(k)=h*(k-1);
end

plot(idok, ertekek(:,1), idok, ertekek(:,2))

%
% Az 1. feladat kódja:

function [idok,ertekek]=elso(h)

% Lépésköz
% h=0.5;

% Végső érték
T=10;

% Lépések száma a [0,T] intervallumon (a kezdetit is beleszámolva)
N = T/h+1;

% Felveszünk egy vektort, amiben tároljuk a lépéseket.
ertekek = zeros(N,1);

% Kezdeti feltétel:
ertekek(1)=1;

% Egy másikban tároljuk az időket (az ábrázoláshoz kell ez majd)
idok = zeros(N,1);

for k=2:N
    ertekek(k) = ertekek(k-1) + h*atan(ertekek(k-1));
    idok(k)=h*(k-1);
end

end