%close all;

%Rugo egyenlet
%x'=y
%y'=-x

%Kezdeti ertekek specialisan:
%x(0)=0
%y(0)=1

%KEP megoldasa:
%x(t) = sin(t)
%y(t) = cos(t)

%VAGY AHOGY A JEGYZETBEN VAN:
%x(t) = cos(t)
%y(t) = -sin(t)

%5 Numerikus modszert vizsgalunk.
%Azt tapasztaljuk, hogy minosegileg maskepp viselkednek a parameterek
%fuggvenyeben.
%Kiszamoljuk a rendszerek energiait es abrazoljuk kulonfele koordinatak
%szerint az egyes megoldasokat.

%Jeloles magyarazat:

%StepSize === h
%NumOfSteps === N (ertsd: 1->N)
%time === [ t_0 t_1 ...t_n... t_N-1 ]
%EEuler === Explicit Euler modszer matrixa
%analytic === Analitikus megoldas matrixa
%IEuler === Implicit Euler modszer matrixa
%SEuler === Semi Euler modszer matrixa
%NRG_*E === A rendszer energiajanak vektora az ido szerint
%NRG_AN === A rendszer energiajanak vektora az ido szerint

%Matrix jeloles:
% Matrix/Vektor elso elemenek indexe a Matlabban 1!!!
%*Euler(a,1) === x === x_n
%*Euler(a,2) === y === y_n
%*Euler(a+1,1) === X === x_n+1
%*Euler(a+1,2) === Y === y_n+1

%FELADAT:
%1) Implementaljuk az egyes modszereket a jeloleseknek megfeleloen
%2) Probaljunk ki mas lepeskoz es lepesszam kombinaciokat! Mit latunk?
StepSize=0.1; %0.1;%0.2%Ez csunya nagy lepeskoz%   0.001;%
NumOfSteps=100; %1000;%50%Ez keves lepes%  100000;%

x_null=0;
y_null=1;
t_null=0;

% Ido kiszamolasa.
%Minden modszer azonos diszkret ido-pontokat hasznal, ezert itt kulon van
%szedve.
time=zeros(NumOfSteps,1);
time(1)=t_null;

for a=1:(NumOfSteps-1)
    time(a+1)=time(a)+StepSize;
end


%Analytic solution
%x(t) = cos(t)
%y(t) = -sin(t)

analytic=zeros(NumOfSteps,2);

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

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

%Gyakorlatilag:
%time=linspace(0,StepSize*NumOfSteps,NumOfSteps) %Valahogy így
%analytic=[sin(time) cos(time)];

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


%Explicit euler method
% (X-x)/h = y
% (Y-y)/h =-x
EEuler=zeros(NumOfSteps,2);

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

for a=1:(NumOfSteps-1)
    EEuler(a+1,1) = EEuler(a,1) + StepSize .* EEuler(a,2);
    EEuler(a+1,2) = EEuler(a,2) + StepSize .* (-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(NumOfSteps,2);

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

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

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


%Semi-Euler method
SEuler = zeros(NumOfSteps,2);

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

for a = 1 : (NumOfSteps-1)
   % (X-x)/h =y
   % (Y-y)/h =-X
    SEuler(a+1,1)  = SEuler(a,1) + StepSize * SEuler(a,2);
    SEuler(a+1,2)  = SEuler(a,2) - StepSize * (SEuler(a,1) + StepSize * SEuler(a,2));
end

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


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

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

figure;
%t-szerint
subplot(2,2,1);
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',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',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',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');