%%%%%%%%%%%%%%%%%%%%%%%%%%%%     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;                                % 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
NRG_AN=analytic(:,1).*analytic(:,1)+analytic(:,2).*analytic(:,2); %energia
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% kezdeti értékek:
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*(b*EEuler(a,2)+EEuler(a,1));

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

%ode45 megoldás:
rugo_ode=@(t,xy) [xy(2); -b*xy(2)-xy(1)];
options=odeset('RelTol', 1e-6, 'AbsTol', 1e-6 );
[timeode45, ode45mg]=ode45(rugo_ode, [t_null, t_max], [x_null, y_null], options);
NRG_ode45=sqrt(ode45mg(:,1).*ode45mg(:,1)+ode45mg(:,2).*ode45mg(:,2));

figure
%t-xszerint
subplot(2,2,1);
%verzió 1:
plot(time,analytic(:,1),'black');
hold on;
plot(time, EEuler(:,1),'r');
hold off;

% verzió 2:
plot(time,analytic(:,1),'black',time, EEuler(:,1),'r', timeode45, ode45mg(:,1), 'g');



xlabel('t');
ylabel('x');
axis equal
legend('Analitikus','Explicit Euler', 'ode45');

subplot(2,2,2);
plot(analytic(:,1), analytic(:,2), 'black',EEuler(:,1), EEuler(:,2),'r', ode45mg(:,1), ode45mg(:,2), 'g');
xlabel('x');
ylabel('y');
axis equal
legend('Analitikus','Explicit Euler', 'ode45');

subplot(2,2,3);
plot(time,NRG_AN,'black',time, NRG_EE,'r', timeode45, NRG_ode45, 'g');
xlabel('t');
ylabel('energia');
axis equal
legend('Analitikus','Explicit Euler', 'ode45');

subplot(2,2,4);
plot3(analytic(:,1),analytic(:,2), time, 'black', EEuler(:,1),EEuler(:,2), time,'r', ode45mg(:,1), ode45mg(:,2), timeode45, 'g');
xlabel('x');
ylabel('y');
zlabel('t');
legend('Analitikus','Explicit Euler', 'ode45');







