%% Metody Rungeho-Kutty 2. stupně
%
% Obecná metoda: y_{i+1}=y_i+h(b_1k_1+b_2k_2), k_1=f(x_i,y_i), k_2=f(x_i+c\,h,y_i+a\,h\,k_1),
% Podmínky řádu: b_1+b_2=1, a=c, a\,b_2=\frac{1}{2}.
% Po dosazení: y_{i+1}=y_i+h[(1-b)f(x_i,y_i)+bf(x_i+ah,y_i+ahf(x_i,y_i))], ab=\frac{1}{2}
% Modifikovaná Eulerova metoda: a=\frac{1}{2},b=1,
% y_{i+1}=y_i+hf(x_i+h,y_i+hf(x_i,y_i))
% Příklad 1: y'=-\frac{y}{x},\ y(1)=1

f=@(x,y)(-y./x);

% Použití modifikované Eulerovy metody:

a=0.5;b=1;      % parametry metody
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
    k1=f(x,y);
    k2=f(x+a*h,y+a*h*k1);
    x1=x+h;
    y1=y+h*((1-b)*k1+b*k2);
    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 modifikovanou Eulerovou 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 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('Modif. E metoda','EE metoda','IE metoda','lich metoda','presn.reseni')

%
% Použití modifikované Eulerovy metody - 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
    k1=f(x,y);
    k2=f(x+a*h,y+a*h*k1);
    x1=x+h;
    y1=y+h*((1-b)*k1+b*k2);
    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 modifikovanou Eulerovou metodou.')

% 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)

%
% Ralstonova metoda 2. řádu: a=\frac{2}{3},b=\frac{3}{4},
% y_{i+1}=y_i+hf(x_i+h,y_i+hf(x_i,y_i))
% Příklad 1: y'=-\frac{y}{x},\ y(1)=1

f=@(x,y)(-y./x);

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

a=2/3;b=3/4;    % parametry metody
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
    k1=f(x,y);
    k2=f(x+a*h,y+a*h*k1);
    x1=x+h;
    y1=y+h*((1-b)*k1+b*k2);
    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 Ralstonovu metodou 2. radu.')
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('Ralstonova metoda','EE metoda','IE metoda','lich metoda','presn.reseni')

%
% Použití Ralstonovy 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
    k1=f(x,y);
    k2=f(x+a*h,y+a*h*k1);
    x1=x+h;
    y1=y+h*((1-b)*k1+b*k2);
    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 Ralstonovu metodou 2. radu.')

% 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)

%
% Heunova metoda: a=1,b=\frac{1}{2},
% y_{i+1}=y_i+hf(x_i+h,y_i+hf(x_i,y_i))
% Příklad 1: y'=-\frac{y}{x},\ y(1)=1

f=@(x,y)(-y./x);

% Použití Heunovy metody:

a=1;b=0.5;      % parametry metody
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
    k1=f(x,y);
    k2=f(x+a*h,y+a*h*k1);
    x1=x+h;
    y1=y+h*((1-b)*k1+b*k2);
    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 Heunovou 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 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('Heunova metoda','EE metoda','IE metoda','lich metoda','presn.reseni')

%
% Použití Heunovy metody - 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
    k1=f(x,y);
    k2=f(x+a*h,y+a*h*k1);
    x1=x+h;
    y1=y+h*((1-b)*k1+b*k2);
    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 Heunovou metodou.')

% 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)