%% Vícekrokové metody
%
%% Adamsovy-Bashforthovy (AB) metody - explicitní
% s+1 kroková metoda: y_{i+1}=y_i+h\sum\limits_{j=0}^sb_j\Delta^jf_{i-j}
% f_k=f(x_k,y_k)=y'_k
% \Delta^0f_k=f_k,\, \Delta f_k=\Delta^1 f_k=f_{k+1}-f_k,\dots,
% \Delta^{r+1} f_k=\Delta^{r}f_{k+1}-\Delta^{1}f_k
% j:      0       1       2       3         4           5
% b_j:    1     1/2   5/12    3/8  251/720  95/288
%
%% Adamsovy-Moultonovy (AM) metody - implicitní
% s+1 kroková metoda: y_{i+1}=y_i+h\sum\limits_{j=-1}^sc_{j+1}\Delta^{j+1}f_{i-j}
% j:      0       1       2        3          4            5
% c_j:    1     -1/2  -1/12  -1/24  -19/720  -31/160
%
%% Obecná vícekroková metoda:
% y_{i+1}=\sum\limits_{j=0}^s a_{j}y_{i-j}+h\sum\limits_{j=-1}^s b_{j}y'_{i-j}
%
%% Formule řádu r:
%% \sum\limits_{j=0}^sa_j=1
%% \sum\limits_{j=0}^s (-j)^ka_j+k\sum\limits_{j=-1}^s(-j)^{k-1}b_j=1,\ k=1,\dots,r
% Metoda (formule) se nazývá konzistentní, jestliže je řádu alespoň 1.
%% Úlohy: určete řád některých AB a AM metod
%
%% Příklady: testování některých AB a AM metod
% Příklad 1:
% dvoukroková AB metoda y_{i+1}=y_i+h[\frac{3}{2}f_i-\frac{1}{2}f_{i-1}]
% rovnice y'=-y^2,\,y(1)=1,

f=@(x,y) -y.^2; % funkce
x0=1; y0=1;     % pocatecni hodnoty
h=0.5;          % krok - zkuste 0.5, 0.1
xn=5;           % koncovy bod
x1=x0+h;
y1=1/x1;        % 2. pocatecni hodnota (presne)
xx=[x0,x1]; 
yy=[y0,y1];     % vysledne priblizne reseni
while x1<xn
    %DOPLNIT
end
plot(xx,yy,xx,1./xx,'--')   % graf priblizneho reseni a presneho reseni
title('Reseni dif. rovnice dvoukrokovou AB metodou.')
legend('num. reseni','presn.reseni')

% Rozdíl mezi přesným a přibližným češením:

R=max(abs(yy-1./xx));
fprintf('Max. chyba reseni je %10.8f\n',R)

% Porovnání s lichoběžníkovou metodou

[xl,yl]=lm(f,1,xn,1,h);
plot(xx,yy,xl,yl,xx,1./xx,'g--')
title('Porovnani s lichobb. metodou.')
legend('AB metoda','lich. metoda','presne reseni')

%
% Příklad 1:
% tříkroková AB metoda y_{i+1}=y_i+\frac{h}{12}[23f_i-16f_{i-1}+5f_{i-2}]
% rovnice y'=-y^2,\,y(1)=1,

f=@(x,y) -y.^2; % funkce
x0=1; y0=1;     % pocatecni hodnoty
h=0.5;          % krok - zkuste 0.5, 0.1
xn=5;           % koncovy bod
x1=x0+h;
y1=1/x1;        % 2. pocatecni hodnota
x2=x1+h;
y2=1/x2;        % 3. pocatecni hodnota
xx=[x0,x1,x2]; 
yy=[y0,y1,y2];  % vysledne priblizne reseni
while x2<xn
    %DOPLNIT
end
plot(xx,yy,xx,1./xx,'--')   % graf priblizneho reseni a presneho reseni
title('Reseni dif. rovnice dvoukrokovou AB metodou.')
legend('num. reseni','presn.reseni')

% Porovnání s lichoběžníkovou metodou a výsledkem funkce ode45

[xl,yl]=lm(f,1,xn,1,h);
[x45,y45]=ode45(f,xx,1);
plot(xx,yy,xl,yl,x45,y45,xx,1./xx,'--')
title('Porovnani s lichobb. metodou.')
legend('AB2 metoda','lich. metoda','ode45','presne reseni')
plot(xx,yy-1./xx,xx,yl-1./xl,x45,y45-1./x45);
title('Rozdil mezi priblizneho a presneho reseni.')
legend('AB2 metoda','lich. metoda','ode45')

% Rozdíl mezi přesným a přibližným češením:

R=max(abs(yy-1./xx));
fprintf('Max. chyba reseni pro AB2 metodu je %10.8f\n',R)
Rl=max(abs(yl-1./xl));
fprintf('Max. chyba reseni pro lich. metodu je %10.8f\n',Rl)
R45=max(abs(y45-1./x45));
fprintf('Max. chyba reseni pro ode45 je %10.8f\n',R45)

%

function [xx,yy]=lm(f,x0,xn,y0,h)
% function [xx,yy]=lm(f,x0,xn,y0,h)
% lichobeznikova  metoda pro rovnici y'=f(x,y)
% f  - anonymni funkce promennych x,y
% x0 - zacatek intervalu
% xn - konec intervalu
% y0 - pocatecni hodnota
% h  - krok
% nelinearni rovnice se resi Newtonovou metodou

presnost=0.000001;
x=x0;
y=y0;
xx=x0;
yy=y0;
fs_y=diff(sym(f),'y');                          % derivace f podle y symbolicky
fy=matlabFunction(fs_y,'Vars',{'x','y'});   % derivace f podle y
while x<xn
    fi=f(x,y);      % hodnota funkce f v predchoyim bode
    z0=y+h*f(x,y);
    yi=y;
    xi=x;
    x=x+h;
    z0=y+h*fi;      % pocatecni aproximace reseni nelin. rovnice je hodnota ziskana s pomoci EEM
    y=z0;
    z1=z0-(z0-yi-0.5*h*(fi+f(x,y)))/(1-0.5*h*fy(x,y));    % prvni aproximace rseni nelin. rovnice
    while abs(z1-z0)>0.001                                  % Newtonova metoda
        y=z1;
        z0=z1;
        z1=z0-(z0-yi-0.5*h*(fi+f(x,y)))/(1-0.5*h*fy(x,y));
    end
    y=z1;
    xx=[xx,x];
    yy=[yy,y];
end
end

%