%% Metoda Taylorova rozvoje
% 1.řádu:    y_{i+1}=y_i+hf(x_i,y_i)
% 2.řádu:    y_{i+1}=y_i+hf(x_i,y_i)₊\frac{1}{2}h^2[f_x(x_i,y_i)+f_y(x_i,y_i)f(x_i,y_i)]
% 3. řádu:    y_{i+1}=y_i+hf(x_i,y_i)₊\frac{1}{2}h^2[f_x(x_i,y_i)+f_y(x_i,y_i)f(x_i,y_i)]
% +\frac{1}{6}h^3[f_{xx}+2f_{xy}f+f_{yy}f^2+f_xf_y+f_y^2f](x_i.y_i)
% Příklad 1: y'=-\frac{y}{x},\ y(1)=1

f=@(x,y)(-y./x);
fs=sym(f);
fs_x=diff(fs,'x');
fs_y=diff(fs,'y');
f_x=matlabFunction(fs_x)
f_y=matlabFunction(fs_y)

% POZOR! Předchozí výsledek může být funkce jen proměnné x.
% Správný postup:

f_x=matlabFunction(fs_x,'Vars',{'x','y'})
f_y=matlabFunction(fs_y,'Vars',{'x','y'})

% Výpočet druhých parciálních derivací:

fs_xx=diff(fs_x,'x');
fs_xy=diff(fs_x,'y');
fs_yy=diff(fs_y,'y');
f_xx=matlabFunction(fs_xx,'Vars',{'x','y'});
f_xy=matlabFunction(fs_xy,'Vars',{'x','y'});
f_yy=matlabFunction(fs_yy,'Vars',{'x','y'});

%
% Použití metody 2. řádu:

x0=1; y0=1;     % pocatecni hodnoty
h=0.5;          % krok
xn=5;           % koncovy bod
xx=x0; yy=y0;   % vysledne priblizne reseni
x=x0;y=y0;
while x<xn
    x1=x+h;
    y1=y+h*f(x,y)+h^2/2*(f_x(x,y)+f_y(x,y)*f(x,y));
    xx=[xx,x1];
    yy=[yy y1];
    x=x1;
    y=y1;
end
plot(xx,yy,xx,1./xx,'--')   % graf priblizneho reseni a presneho reseni
title('Reseni dif. rovnice metodou Taylorova rozvoje.')
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 Eulerovými metodami

[x1,y1]=ee(f,x0,xn,y0,h);
[x2,y2]=ie(f,x0,xn,y0,h);
[x3,y3]=lm(f,x0,xn,y0,h);
plot(xx,yy,x1,y1,x2,y2,x3,y3,xx,1./xx,'g--')
title('Porovnani s Eulerovymi metodami.')
legend('Taylor. rozvoj','EE metoda','IE metoda','lich metoda','presn.reseni')

%
% Použití metody 2. řádu - menší krok:

x0=1; y0=1;     % pocatecni hodnoty
h=0.1;          % krok
xn=5;           % koncovy bod
xx=x0; yy=y0;   % vysledne priblizne reseni
x=x0;y=y0;
while x<xn
    x1=x+h;
    y1=y+h*f(x,y)+h^2/2*(f_x(x,y)+f_y(x,y)*f(x,y));
    xx=[xx,x1];
    yy=[yy y1];
    x=x1;
    y=y1;
end
plot(xx,yy,xx,1./xx,'--')   % graf priblizneho reseni a presneho reseni
legend('num. reseni','presn.reseni')
title('Resení dif. rovnice metodou Taylorova rozvoje.')

% 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 Eulerovými metodami - jen hodnoty chyb

[x1,y1]=ee(f,x0,xn,y0,h);
[x2,y2]=ie(f,x0,xn,y0,h);
[x3,y3]=lm(f,x0,xn,y0,h);
R1=max(abs(y1-1./x1));
R2=max(abs(y2-1./x2));
R3=max(abs(y3-1./x3));
fprintf('Max. chyba reseni pro EE metodu je %10.8f\n',R1)
fprintf('Max. chyba reseni pro IE metodu je %10.8f\n',R2)
fprintf('Max. chyba reseni pro lich. metodu je %10.8f\n',R3)

%
% Použití metody 3. řádu:

x0=1; y0=1;     % pocatecni hodnoty
h=0.5;          % krok
xn=5;           % koncovy bod
xx=x0; yy=y0;   % vysledne priblizne reseni
x=x0;y=y0;
while x<xn
    x1=x+h;
    y1=y+h*f(x,y)+h^2/2*(f_x(x,y)+f_y(x,y)*f(x,y));
    y1=y1+h^3/6*(f_xx(x,y)+2*f_xy(x,y)*f(x,y)+f_yy(x,y)*f(x,y)^2+f_x(x,y)*f_y(x,y)+f_y(x,y)^2*f(x,y));
    xx=[xx,x1];
    yy=[yy y1];
    x=x1;
    y=y1;
end
plot(xx,yy,xx,1./xx,'--')   % graf priblizneho reseni a presneho reseni
title('Reseni dif. rovnice metodou Taylorova rozvoje.')
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 Eulerovými metodami

[x1,y1]=ee(f,x0,xn,y0,h);
[x2,y2]=ie(f,x0,xn,y0,h);
[x3,y3]=lm(f,x0,xn,y0,h);
plot(xx,yy,x1,y1,x2,y2,x3,y3,xx,1./xx,'g--')
title('Porovnani s Eulerovymi metodami.')
legend('Taylor. rozvoj','EE metoda','IE metoda','lich metoda','presn.reseni')

%
% Použití metody 3. řádu - menší krok:

x0=1; y0=1;     % pocatecni hodnoty
h=0.1;          % krok
xn=5;           % koncovy bod
xx=x0; yy=y0;   % vysledne priblizne reseni
x=x0;y=y0;
while x<xn
    x1=x+h;
    y1=y+h*f(x,y)+h^2/2*(f_x(x,y)+f_y(x,y)*f(x,y));
    y1=y1+h^3/6*(f_xx(x,y)+2*f_xy(x,y)*f(x,y)+f_yy(x,y)*f(x,y)^2+f_x(x,y)*f_y(x,y)+f_y(x,y)^2*f(x,y));
    xx=[xx,x1];
    yy=[yy y1];
    x=x1;
    y=y1;
end
plot(xx,yy,xx,1./xx,'--')   % graf priblizneho reseni a presneho reseni
legend('num. reseni','presn.reseni')
title('Resení dif. rovnice metodou Taylorova rozvoje.')

% 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 Eulerovými metodami - jen hodnoty chyb

[x1,y1]=ee(f,x0,xn,y0,h);
[x2,y2]=ie(f,x0,xn,y0,h);
[x3,y3]=lm(f,x0,xn,y0,h);
R1=max(abs(y1-1./x1));
R2=max(abs(y2-1./x2));
R3=max(abs(y3-1./x3));
fprintf('Max. chyba reseni pro EE metodu je %10.8f\n',R1)
fprintf('Max. chyba reseni pro IE metodu je %10.8f\n',R2)
fprintf('Max. chyba reseni pro lich. metodu je %10.8f\n',R3)