%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         4.gyakorlat         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nemlin. din. rsz-ek alapjai %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        2014.10.08/09        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     Numerikus módszerek     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Paraméterek:
h=0.1;                      % lepeskoz
N=100;                      % lepesszam

x_null=0;                   % kezdeti feltételek
y_null=1;
t_null=0;

% Idõ kiszámolása.
% Minden módszer azonos diszkrét idõpontokat használ, ezért itt külon van
% szedve.
time=zeros(N,1);            % N*1-es mátrix végig 0 elemekkel
time(1)=t_null;             % kezdeti érték megadása
for a=1:(N-1)               % "idõ mátrix" feltöltése
    time(a+1)=time(a)+h;
end
% FELADAT
% Adott a következõ egyenlet: x"+x = 0
% Átírva 2 elsõfokúra:
% x'= y
% y'=-x adódik.

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


% ode45:
rugo_ode=@(t,xy) [xy(2); -xy(1)];
[timeode45, ode45mg]=ode45(rugo_ode, [t_null, h*N], [x_null, y_null]);

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


% Analitikus megoldás:

analytic=zeros(N,2);
analytic(1,1)=x_null;
analytic(1,2)=y_null;

for a=1:(N-1)
    analytic(a+1,1) = sin(time(a+1));
    analytic(a+1,2) = cos(time(a+1));
end

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



% Explicit Euler módszer
% (X-x)/h = y
% (Y-y)/h =-x

EEuler=zeros(N,2);

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

for a=1:(N-1)
    EEuler(a+1,1) = EEuler(a,1) + h .* EEuler(a,2);
    EEuler(a+1,2) = EEuler(a,2) + h .* (-EEuler(a,1));
end

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



%Implicit Euler method
% (X-x)/h = Y
% (Y-y)/h =-X
IEuler = zeros(N,2);

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

for a = 1 : (N-1)
    IEuler(a+1,1)  = (IEuler(a,1)  + (h *IEuler(a,2))) / (1+h*h);
    IEuler(a+1,2)  = (IEuler(a,2)- (h * IEuler(a,1))) / (1+h*h);
end

% Euklideszi norma:
NRG_IE=IEuler(:,1).*IEuler(:,1)+IEuler(:,2).*IEuler(:,2);



% Szemi-implicit Euler:
   % (X-x)/h =y
   % (Y-y)/h =-X
SEuler = zeros(N,2);

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

for a = 1 : (N-1)
    SEuler(a+1,1)  = SEuler(a,1) + h * SEuler(a,2);
    SEuler(a+1,2)  = SEuler(a,2) - h * (SEuler(a,1) + h * SEuler(a,2));
end

% Euklideszi norma:
NRG_SE=SEuler(:,1).*SEuler(:,1)+SEuler(:,2).*SEuler(:,2);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        Kirajzoltatás        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% az implicit és szemi-implicit Euler ki van comment-elve

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',time,IEuler(:,1),'green',time,SEuler(:,1),'blue', timeode45, ode45mg(:,1),'magenta','DisplayName', 't szerint');
xlabel('t');
ylabel('x');
axis equal
legend('Anal', 'EE', 'IE', 'SE', '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',IEuler(:,1),IEuler(:,2),'green',SEuler(:,1),SEuler(:,2),'blue', ode45mg(:,1), ode45mg(:,2), 'magenta','DisplayName', 'x-y szerint')
xlabel('x');
ylabel('y');
axis equal
legend('Anal', 'EE', 'IE', 'SE', 'ode45');

%Energiak
subplot(2,2,3);
%plot(time,NRG_AN,'black',time,NRG_EE,'red');
plot(time,NRG_AN,'black',time,NRG_EE,'red',time,NRG_IE,'green',time,NRG_SE,'blue', timeode45, NRG_ode45,'magenta','DisplayName', 'Energiak')
xlabel('t');
ylabel('Energia');
axis equal
legend('Anal', 'EE', 'IE', 'SE', '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',time,IEuler(:,1),IEuler(:,2),'green',time,SEuler(:,1),SEuler(:,2),'blue', timeode45, ode45mg(:,1), ode45mg(:,2),'magenta','DisplayName', 't-x-y szerint')
grid on
xlabel('t');
ylabel('x');
zlabel('y');
legend('Anal', 'EE', 'IE', 'SE', 'ode45');