% Demonstrace Eulerovy metody a Runge-Kuttovy metody ctvrteho radu pro reseni obycejnych diferencialnich rovnic
%
% Budeme resit harmonicky oscilator, tedy rovnici y''+omega^2*y=0
% Napred musime rovnici druheho radu prevest na soustavu rovnic prvniho radu
% takto
%   y'=z
%   z'=-omega^2*y
% To jeste prepiseme jako ze u'=f(t,u), tedy u=[y,z]
% a dale budeme chtit pocatecni stav takovy, ze rychlost na pocatku je v0=0 a
% vychylka y0 je nenulova

v0=0; % pocatecni rychlost je nula, tedy y'(0)=v(0)=v0=z1(0)=0
y0=1; % pocatecni vychylka je nenulova
T=20; % budeme integrovat 20 sekund
h=0.5; % pouzijeme krok h
%h=0.1;
%h=0.01;
omega=0.5; % nastavime uhlovou rychlost

u1=[y0,v0]; % nastavime na zacatek promennou u
u2=[y0,v0]; 

f=inline(sprintf('[u(2),-%g^2*u(1)]',omega),'t','u'); % definice funkce f, plati ze u'=f(t,u)


figure;  % otevreme novy obrazek, nastavime osy a popisky
subplot(211);
title('Poloha');
xlabel('t [s]');
ylabel('y [m]');
axis([0 T -1.5 1.5]);
subplot(212);
title('Rychlost');
xlabel('t [s]');
ylabel('v [m/s]');
axis([0 T -1 1]);

% presne reseni: jen pro kresleni
xx=(0:T/100:T);
yex = cos(omega*xx);
dyex=-omega*sin(omega*xx);
subplot(211);   % nakreslime presne reseni
hold on;
plot(xx,yex,'k');
hold off;
subplot(212);
hold on;
plot(xx,dyex,'k');
hold off;

pocetKroku=T/h; % celkovy pocet kroku
i=0;
while 1,
   if i>=pocetKroku, % skoncime kdyz jsme dosahli potrebneho mnozstvi kroku
       break;
   end;
   
   t=i/pocetKroku*T; % spocteme cas odpovidajici kroku
   
   subplot(211);
   hold on;% nastavime grafy aby se nemazaly, ale aby se jen pridavaly objekty do existujicich grafu
   plot(t,u1(1),'b.'); % zobrazime aktualni polohu
   plot(t,u2(1),'r+'); % zobrazime aktualni polohu
   hold off;
   subplot(212);
   hold on;
   plot(t,u1(2),'b.'); % zobrazime aktualni rychlost
   plot(t,u2(2),'r+'); % zobrazime aktualni rychlost
   hold off;
   
   % spocteme novou hodnotu promenne u
 
   % Euler
   u1=u1+h*f(t,u1); % spocteme novou hodnotu promenne u
   
   % Runge-Kutta 4
   k1 = f(t,u2);
   k2 = f(t+h/2, u2+h/2*k1);
   k3 = f(t+h/2, u2+h/2*k2);
   k4 = f(t+h, u2+h*k3);
   u2=u2+h/6*(k1+2*k2+2*k3+k4);
   
   i=i+1; % zvetsime krok o jednicku

end;

