%%%%%%%%%%%%%%%%%%%%%%%%%%%%     3. gyakorlat    %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%      2015.02.26.    %%%%%%%%%%%%%%%%%%%%%%%%%%%

% FELADAT:
% Adott a következõ egyenlet: x"+x = 0
% Átírva 2 elsõfokúra:
% x'= y
% y'=-x adódik.

% csillapított egyenlet: x"+b*x'+x=0
% x'=y
% y'=-b*x'-x=-b*y-x

% A DER analitikus megoldása:
% dx(t) = cos(t)
% dy(t) = -sin(t)

% 0. Paraméterek 

t_max=100;                          % idõ vége
t_null=0;                           % idõ eleje
h=0.1;                              % lépésköz
time=t_null:h:t_max;                % idõpillanatok

x_null=0;                           % kezdeti feltételek megadása
y_null=1;
b=0.5;                                % csillapítás

% I. Analitikus megoldás:

analytic=zeros(length(time),2);     % time hosszú*2 0-kból álló mátrix
analytic(1,1)=x_null;               % x kezdeti érték beállítása a mátrixban
analytic(1,2)=y_null;               % y kezdeti érték beállítása a mátrixban

for a=1:length(time)                % analitikus megoldás, diszkrét idõ
    analytic(a,1) = sin(time(a));
    analytic(a,2) = cos(time(a));
end

% Euklideszi norma:

NRG_AN=analytic(:,1).*analytic(:,1)+analytic(:,2).*analytic(:,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% II. Numerikus módszer: Explicit Euler

% (X-x)/h = y
% (Y-y)/h =-x


EEuler=zeros(length(time),2);

EEuler(1,1)=x_null;
EEuler(1,2)=y_null;

for a=1:length(time)-1
    EEuler(a+1,1) = EEuler(a,1) + h .* EEuler(a,2);
%   EEuler(a+1,2) = EEuler(a,2) + h .* (-EEuler(a,1));
    EEuler(a+1,2) = EEuler(a,2) + h .* (-b*EEuler(a,2)-EEuler(a,1)); %csillapított    
end

% Euklideszi norma:

NRG_EE=EEuler(:,1).*EEuler(:,1)+EEuler(:,2).*EEuler(:,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ode45

% rugo_ode=@(t,xy) [xy(2); -xy(1)];
 rugo_ode=@(t, xy) [xy(2); -b*xy(2)-xy(1)]; %csillapított
% [timeode45, ode45mg]=ode45(rugo_ode, [t_null, t_max], [x_null, y_null]);
 %options=odeset('RelTol', 1e-6, 'AbsTol', 1e-9 ); %az ode45 pontosításására (alap: 1e-3, 1e-6) %inkább az AbsTol számít, ha valamihez kell közel lenni (egyens pont), és a RelTol ha egymáshoz (periódus)?
% [timeode45, ode45mg]=ode45(rugo_ode, [t_null, t_max], [x_null, y_null], options);

% másik módszer: függvény meghívás
% [timeode45 ode45mg] = rugo_numgyakhoz(t_null, t_max, x_null, y_null,b);
 [timeode45, ode45mg] = ode45(@(timeode45, ode45mg) rugo_egyenlet(timeode45, ode45mg, b),[t_null, t_max], [x_null, y_null]);

% Euklideszi norma:
NRG_ode45=ode45mg(:,1).*ode45mg(:,1)+ode45mg(:,2).*ode45mg(:,2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Kirajzoltatás   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;
%t-szerint
subplot(2,2,1);
%plot(time, analytic(:,1),'black',time,EEuler(:,1),'red');
plot(time,analytic(:,1),'black',time,EEuler(:,1),'red', timeode45, ode45mg(:,1),'green','DisplayName', 't szerint');
xlabel('t');
ylabel('x');
axis equal
legend('Anal', 'EE', 'ode45');

%x-y szerint
subplot(2,2,2)
%plot(analytic(:,1),analytic(:,2),'black',EEuler(:,1),EEuler(:,2),'red');
plot(analytic(:,1),analytic(:,2),'black',EEuler(:,1),EEuler(:,2),'red', ode45mg(:,1), ode45mg(:,2), 'green','DisplayName', 'x-y szerint')
xlabel('x');
ylabel('y');
axis equal
legend('Anal', 'EE', 'ode45');

%Energiak
subplot(2,2,3);
%plot(time,NRG_AN,'black',time,NRG_EE,'red');
plot(time,NRG_AN,'black',time,NRG_EE,'red', timeode45, NRG_ode45,'green','DisplayName', 'Energiak')
xlabel('t');
ylabel('Energia');
axis equal
legend('Anal', 'EE', 'ode45');

%t-x-y szerint
subplot(2,2,4);
%plot3(time,analytic(:,1),analytic(:,2),'black',time,EEuler(:,1),EEuler(:,2),'red');
plot3(time,analytic(:,1),analytic(:,2),'black',time,EEuler(:,1),EEuler(:,2),'red', timeode45, ode45mg(:,1), ode45mg(:,2),'green','DisplayName', 't-x-y szerint')
grid on
xlabel('t');
ylabel('x');
zlabel('y');
legend('Anal', 'EE', 'ode45');

